Confusion over density-density correlation function I have a question regarding what may perhaps be quite simple definitions in classical statistical physics. If one considers a set of $n$ random variables with an $n$-dimensional probability density function $p(u_1,\ldots,u_n)$, then the moments of order $k = \sum_{i=1}^n k_i$ are defined via
$\langle u_1^{k_1}\ldots u_n^{k_n} \rangle = \int_{-\infty}^\infty\ldots\int_{-\infty}^\infty u_1^{k_1}\ldots u_n^{k}\,p(u_1,\ldots,u_n)\,du_1\ldots du_n$.
Further, expressions of the form $\langle u(\mathbf{x}_1)u(\mathbf{x}_2) \rangle$ are known as correlation functions. This all seems simple enough, until I came across a concrete example. Solid crystals are in some sense defined by the fact that their density-density correlation functions satisfy the following property:
$\lim_{|\mathbf{x} - \mathbf{x'}|}\langle \rho(\mathbf{x})\rho(\mathbf{x}') \rangle = f(\mathbf{x}-\mathbf{x}')$, where $f$ is periodic in the chosen set of basis vectors. However, based on the above definitions, I propose that
$\langle \rho(\mathbf{x})\rho(\mathbf{x}') \rangle \stackrel{?}{=} \int\int\rho(\mathbf{x})\rho(\mathbf{x}')p(\mathbf{x},\mathbf{x}')'\,d^3\mathbf{x}\,d^3\mathbf{x}'$,
where the integrals are taken over all space.
My questions are therefore


*

*Is my understanding of the density-density correlation function correct? 

*If not, what should I understand by $\langle \rho(\mathbf{x})\rho(\mathbf{x}') \rangle$?

*If so, how can the resultant integral be a function of space; and

*how does one understand $p$ in this setting?


I'm sorry if this is more straightforward than I am making it, but I am quite new to statistical mechanics and struggling to find my bearings.
 A: Density function (normalied) of a $N$-particle classical system with cordinates $\{\vec{r}_{i}^{}\}$ can be defined as :
$$\rho(\vec{r})=\frac{1}{N}\sum_{i=1}^{N}\delta_{}^{(3)}(\vec{r}-\vec{r}_{i}^{}).$$
And the density-density correlation function (in canonical ensemble) is defined as :
$$\langle\rho(\vec{r}) \rho(\vec{r}')\rangle = \frac{\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}\rho(\vec{r}) \rho(\vec{r}')e^{-\beta\mathcal{H}(\{\vec{r}_{i}^{}\},\{\vec{p}_{i}^{}\})}}{\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}e^{-\beta\mathcal{H}(\{\vec{r}_{i}^{}\},\{\vec{p}_{i}^{}\})}}.$$
$\textbf{Addendum (Field theoretic generating function technology):}$
To recast the problem in field theoretic terms (to obtain a statistical field theory) for simplicity consider :
$$\mathcal{H}(\{\vec{r}_{i}^{}\},\{\vec{p}_{i}^{}\})=\sum_{i=1}^{N}\frac{\vec{p}_{i}^{}\cdot\vec{p}_{i}^{}}{2m}+\frac{1}{2}\sum_{i\neq j=1}^{N}\mathcal{V}(\vec{r}_{i}^{},\vec{r}_{j}^{}).$$ 
Now define a density moment generating function as :
$$\mathcal{Z}[\{\phi(\vec{r})\}] = \frac{\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}e^{i\int d^{3}_{}\vec{r}\phi(\vec{r})\rho(\vec{r})}_{}e^{-\beta\mathcal{H}(\{\vec{r}_{i}^{}\},\{\vec{p}_{i}^{}\})}}{\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}e^{-\beta\mathcal{H}(\{\vec{r}_{i}^{}\},\{\vec{p}_{i}^{}\})}}.$$
Using which we can calculated density-density correlation functions to all order as :
$$\langle\rho(\vec{r}_{1}^{})\cdots\
\rho(\vec{r}_{n}^{})\rangle=(-i)^{n}_{}\frac{\delta}{\delta \phi(\vec{r}_{1}^{})}\cdots\frac{\delta}{\delta \phi(\vec{r}_{n}^{})}\mathcal{Z}[\{\phi(\vec{r})\}]\Big|_{\{\phi(\vec{r})=0\}}^{}.$$
Recast the moment generating function as (self interaction of particles need to be regularized somehow (denoted as $\Delta$ here)) :
$$\mathcal{Z}[\{\phi(\vec{r})\}] = \frac{\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}e^{i\int d^{3}_{}\vec{r}\phi(\vec{r})\rho(\vec{r})}_{}e^{-\beta\big[\sum_{i=1}^{N}\frac{\vec{p}_{i}^{}\cdot\vec{p}_{i}^{}}{2m}+\frac{N^{2}_{}}{2}\int\int d^{3}_{}\vec{r} d^{3}_{}\vec{r}'\rho(\vec{r})\mathcal{V}(\vec{r},\vec{r}')\rho(\vec{r}')+\Delta[\rho(\vec{r})]\big]}}{\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}e^{-\beta\mathcal{H}(\{\vec{r}_{i}^{}\},\{\vec{p}_{i}^{}\})}}.$$
Now we can introduce the functional delta identity : 
$$\int\mathcal{D}[\sigma(\vec{r})]\delta[\sigma(\vec{r})-\rho(\vec{r})]=1$$
inside the generating function to get :
$$\mathcal{Z}[\{\phi(\vec{r})\}] = \frac{\int\mathcal{D}[\sigma(\vec{r})]e^{i\int d^{3}_{}\vec{r}\phi(\vec{r})\sigma(\vec{r})}_{}e^{-\beta\big[\frac{N^{2}_{}}{2}\int\int d^{3}_{}\vec{r} d^{3}_{}\vec{r}'\sigma(\vec{r})\mathcal{V}(\vec{r},\vec{r}')\sigma(\vec{r}')+\Delta[\sigma(\vec{r})]\big]}e^{-\beta\Phi[\sigma(\vec{r})]}}{\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}e^{-\beta\mathcal{H}(\{\vec{r}_{i}^{}\},\{\vec{p}_{i}^{}\})}}.$$
Where :
$$e^{-\beta\Phi[\sigma(\vec{r})]}=\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}\delta[\sigma(\vec{r})-\rho(\vec{r})]e^{-\beta\big[\sum_{i=1}^{N}\frac{\vec{p}_{i}^{}\cdot\vec{p}_{i}^{}}{2m}\big]}$$
which can be rewritten using functional Fourier representation of functional delta functional as :
$$e^{-\beta\Phi[\sigma(\vec{r})]}=\int\mathcal{D}[\omega(\vec{r})]\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}e^{-i\int d_{}^{3}\vec{r}\omega(\vec{r})[\sigma(\vec{r})-\rho(\vec{r})]}e^{-\beta\big[\sum_{i=1}^{N}\frac{\vec{p}_{i}^{}\cdot\vec{p}_{i}^{}}{2m}\big]}.$$
Finally everything can be recast as :
$$\mathcal{Z}[\{\phi(\vec{r})\}] = \frac{\int\mathcal{D}[\omega(\vec{r})]\int\mathcal{D}[\sigma(\vec{r})]e^{i\int d^{3}_{}\vec{r}[\phi(\vec{r})-\omega(\vec{r})]\sigma(\vec{r})}_{}e^{-\beta\big[\frac{N^{2}_{}}{2}\int\int d^{3}_{}\vec{r} d^{3}_{}\vec{r}'\sigma(\vec{r})\mathcal{V}(\vec{r},\vec{r}')\sigma(\vec{r}')+\Delta[\sigma(\vec{r})]\big]}e^{-\beta F[\omega(\vec{r})]}}{\mathcal{Z}[\{0\}]}.$$
With : 
$$e^{-\beta F[\omega(\vec{r})]}=\int\dots\int\prod_{i=1}^{N}d_{}^{3}\vec{r}_{i}^{}d_{}^{3}\vec{p}_{i}^{}e^{-\beta\sum_{i=1}^{N}\big[\frac{\vec{p}_{i}^{}\cdot\vec{p}_{i}^{}}{2m}+i\beta_{}^{-1}\omega(\vec{r}_{i}^{})\big]},$$
is the partition function of non-interacting particles in a external imaginary potential.
$\textbf{Probability functional for densities:}$ Finally, on inverse functional Fourier transforming $\mathcal{Z}[\{\phi(\vec{r})\}]$ we get probability functional for densities as :
$$P[\{\sigma(\vec{r})\}]=\frac{e^{-\beta\big[\frac{N^{2}_{}}{2}\int\int d^{3}_{}\vec{r} d^{3}_{}\vec{r}'\sigma(\vec{r})\mathcal{V}(\vec{r},\vec{r}')\sigma(\vec{r}')+\Delta[\sigma(\vec{r})]+\mathbb{V}_{}^{}[\sigma(\vec{r})]\big]}}{\int\mathcal{D}[\sigma(\vec{r})]e^{-\beta\big[\frac{N^{2}_{}}{2}\int\int d^{3}_{}\vec{r} d^{3}_{}\vec{r}'\sigma(\vec{r})\mathcal{V}(\vec{r},\vec{r}')\sigma(\vec{r}')+\Delta[\sigma(\vec{r})]+\mathbb{V}_{}^{}[\sigma(\vec{r})]\big]}}$$
with
$$e^{-\beta\mathbb{V}_{}^{}[\sigma(\vec{r})]}=\int\mathcal{D}[\omega(\vec{r})]e^{-i\int d^{3}_{}\vec{r}\omega(\vec{r})\sigma(\vec{r})}_{}e^{-\beta F[\omega(\vec{r})]}.$$
A: You need to pay attention to exactly what the is the random variable you are averaging over. When we measure the density-density correlation function, we choose which points to measure at, so $\mathbf{x}$ is not a random variable, it is a parameter. The randomly fluctuating quantities are $\rho(\mathbf{x})$, that is we have an (infinite) family of correlated random variables, $\rho$ which we parameterise with $\mathbf{x}$. 
This means that when we calculate averages we must average over all values of $\rho$, and since all the $\rho$s are generally correlated we must do this for every $\rho$. This means, rather than the integral you wrote down, we must perform a functional integral over all possible functions $\rho(\mathbf{x})$
$$
\langle \rho(\mathbf{x})\rho(\mathbf{x}')\rangle = \int \mathcal{D}\rho\; \rho(\mathbf{x})\rho(\mathbf{x}') p[\rho]
$$
where the probability $p$ is now a functional of $\rho$
In practice we normally write $p$ in terms of some action $p = \exp(-S[\rho])$ and normally write the entire integral in terms of some underlying field for which we have a microscopic model, such that $\rho = \psi^\dagger \psi$.
