Van der Vaart states in his book "Asymptotic Statistics" that an empirical process (of the sample $X_1, X_2, \dots, X_n$) is defined by $$G_n = \sqrt n( P_n - P),$$ where $P$ is the true probability distribution in question and $P_n = n^{-1}\sum_{j=1}^n \delta_{X_j}$ with $\delta_x$ denoting the Dirac measure at $x$. This process $G_n$ can be indexed by a set of measurable functions $\mathcal F$ by letting $$G_nf = \sqrt n( P_nf - Pf)$$ for every $f\in\mathcal F$. Here $P_nf := n^{-1}\sum_{j=1}^n f(X_j)$ and $Pf := \int f\,\mathrm dP$.
Such a class of functions is called $P$-Donsker iff the process $G_n$ indexed by $\mathcal F$ converges in $\ell^\infty(\mathcal F)$ to a tight random element $\Psi$. Since $\Psi$ is tight, it is completely determined by its finite-dimensional distributions, i.e. by the distribution of $(\Psi(f_1), \Psi(f_2), \dots, \Psi(f_k))$ for each $k\in\mathbb N$.
Now suppose that $X_1, X_2,\dots, X_n$ are an iid sample from a probability distribution on $\mathbb R^p$. Let $F$ denote the vector of marginal distribution functions, i.e. $F:\mathbb R^p\rightarrow\mathbb R^p$ with $x\mapsto F(x)$, where $F = (F_1, F_2, \dots, F_p)$, and $F_i = n^{-1}\sum_{j=1}^n I_{(-\infty; \langle e_i,X_j\rangle]}$. Here $e_i$ denotes the $i$th canonical basis vector of $\mathbb R^p$ and $I_A$ is the indicator function of the set $A$. So let's put $\mathcal F = \{x\mapsto I_{(-\infty; \langle e_i,x\rangle]} : i = 1,2,\dots, p\}$. That is, $\mathcal F$ is the set of indicator functions. There is a proof in van der Vaart's book that $\mathcal F$ is $P$-Donsker, and the random element $\Psi$ is given by a Brownian bridge.
I tried to figure out the variance-covariance function of $\Psi$ but I don't understand how the covariances are obtained. Let's take $p = 1$ for a second. Then it's well-known that the variance-covariance function of $\Psi$ is given by $$(t,s)\mapsto \begin{bmatrix} F_1(t)(1-F_1(t)) & F_1(s \wedge t)(1- F_1(s\vee t)) \\ F_1(s \wedge t)(1- F_1(s\vee t)) & F_1(s)(1-F_1(s))\end{bmatrix}.$$ But where do the off-diagonal elements come from? In van der Vaart book it is stated that the variance-covariance is given by $Pf_if_j - Pf_iPf_j$. If $p = 1$ there should be no covariances. So I suspect that this result indexes the empirical process not by a function but by the domain $T$ of the random variable, and $s,t\in T$. So if the empirical process is indexed by functions, the correlations in a dimension with different points $s$ and $t$ in the domain are not present in the variance-covariance matrix? This is confusing me very much. Furthermore, in the case $p>1$, the off-diagonal elements (according to the formula from van der Vaart's book) is $F_{ij} - F_iF_j$, where $F_{ij}$ denotes the 2-dimensional marginal distribution function of the $i$th and $j$th component. Does that mean the correlation is a function?
As mentioned in vdV in the paragraph preceding equation 19.2: Applying the central limit theorem, we know that for each $t \in \mathbb R^k$, the random vector whose coordinates are $\sqrt n (F_n(t_i) - F(t_i))$ converges in distribution to a zero mean gaussian with covariance matrix whose $(i, j)^{th}$ coordinate is given by:
$$\text{Cov}\left[\mathbb 1(X \leq t_i), \mathbb 1(X \leq t_j) \right]$$
Computing this covariance we get:
$$ = \mathbb E \left [ \left ( \mathbb 1(X \leq t_i)- F(X \leq t_i)\right ) \left ( \mathbb 1(X \leq t_j)- F(X \leq t_j)\right ) \right ] $$
$$ = \mathbb E \left [ \mathbb 1(X \leq t_i)\mathbb 1(X \leq t_j)\right ]- F(X \leq t_i)F(X \leq t_j) $$
Now can you see where the max comes from?