Relation between the Glauber correlation functions and statistical correlation

As stated in Wikipedia's page:

Correlation or dependence is any statistical relationship, whether causal or not, between two random variables.

Now in "The Quantum Theory of Optical Coherence" Glauber introduces a set of functions he says quantifies correlations:

The field average (3.3) which determines the counting rate of an ideal photodetector is a particular form of a more general type of expression whose properties are of considerable interest. In the more general expression, the fields $$E^{(-)}$$ and $$E^{(+)}$$ are evaluated at different spacetime points. Statistical averages of the latter type furnish a measure of the correlations of the complex fields at separated positions and times. We shall define such a correlation function, $$G^{(1)}$$, for the $$\mathbf{e}$$ components of the complex fields as $$G^{(1)}(\mathbf{r},t;\mathbf{r}',t')=\operatorname{tr}\left\{\rho E^{(-)}(\mathbf{r},t)E^{(+)}(\mathbf{r}',t')\right\}\tag{3.6}$$

He then more generaly defines the $$n$$-th order correlation functions as a function of $$2n$$ spacetime points $$x_1,\dots, x_{2n}$$:

$$G^{(n)}(x_1,\dots, x_n;x_{n+1},\dots, x_{2n})=\operatorname{tr}\left\{\rho E^{(-)}(x_1)\cdots E^{(-)}(x_n)E^{(+)}(x_{n+1})\cdots E^{(+)}(x_{2n})\right\}\tag{3.8}$$

In all the above, $$\rho$$ is the quantum state of the electromagnetic field.

My question here is the following: why these $$G^{(n)}$$ are called correlation functions? How they are connected with the statistical idea of correlation?

I imagine in some sense they should quantify some sort of statistical relationship between measurements associated to photons at each of the spacetime points (detection of the photons, perhaps?), but this is an impression because of the name. I want to understand really why these $$G^{(n)}$$ functions quantify correlations and what are the correlations they are quantifying.

The $$G^{(n)}$$ are the expectation values of the products of $$2n$$ of the $$E(x_i)$$ variables, since $$\mathrm{Tr}(\rho X)$$ is how one writes the expectation value for $$X$$ in a mixed state with density matrix $$\rho$$. In particular $$G^{(1)}$$ is the expectation value of a product of two random variables. Since we can usually set the expectation value of the fields to zero by adding constants, this is the same as the covariance of the two random variables. The covariance is "unnormalized correlation", the common Pearson correlation coefficient is merely the covariance divided by the standard deviations of the two variables.
So the "correlation functions" really are straightforward measures of (unnormalized) correlation coefficients between the random variables $$E^{\pm}(x_i)$$ in the standard sense of statistics.
• Thanks @ACuriousMind. There's something I think I'm still missing. I know that given two random variables one may quantify linear correlations with the Pearson correlation coefficient. In the setting of the question this would be $G^{(1)}(x,y)=\langle E^+(x)E^-(y)\rangle$. But the $G^{(n)}(x_1,\dots x_n,x_{n+1},\dots, x_{2n})$ involves $2n$ random variables. What is the interpretation then? Isn't correlation something associated to just a pair of random variables? Or it is some generalization of the Pearson approach to take into account higher order correlations?