# Analytical Formulation

The $\min$ function appearing in the average distances is approximated in term of Gaussian functions:

$\min_{j\neq i} r_{ij} = \lim_{\gamma \rightarrow +\infty} \sqrt{ -\frac{1}{\gamma} \log_2 f(\gamma;r_{ij})}$

with

$f(\gamma;r_{ij}) =  \sum_{j \neq i} e^{ -\gamma r_{ij}^2 }$

Now, our basic approximation consists in replacing, for $\gamma$ large, the integrals

$\left\langle \Psi \left| \sum_{i=1}^N \delta (\mathbf{r}-\mathbf{r}_i) \left( \sqrt{ -\frac{1}{\gamma} \log_2 f(\gamma;r_{ij})} \right) \right| \Psi \right\rangle \qquad (1)$

by

$\sqrt{ -\frac{1}{\gamma} \log_2 \left\langle \Psi \left| \sum_{i=1}^N \delta(\mathbf{r}-\mathbf{r}_i) f(\gamma;r_{ij}) \right| \Psi \right\rangle } \qquad (2)$

The expectation values of the minimum distances are now given by:

$d_{\sigma \sigma}(\mathbf{r}) \sim_{\{\gamma {\rm ~large}\}} \sqrt{-\frac{1}{\gamma} \log_2 {\bar f}_{\sigma \sigma}(\gamma;\mathbf{r})}$

$d_{\sigma \bar{\sigma}} (\mathbf{r}) \sim_{\{\gamma {\rm ~large}\}} \sqrt{-\frac{1}{\gamma} \log_2 {\bar f}_{\sigma \bar{\sigma}} (\gamma;\mathbf{r})}$

with the two-electron integrals:

${\bar f}_{\sigma \sigma}(\gamma;\mathbf{r}) = \left \langle \Psi \left| \sum_{i=1}^{N} \delta(\mathbf{r}-\mathbf{r}_i) \sum_{j \neq i ; \sigma_i = \sigma_j}^{N} e^{ -\gamma |\mathbf{r}_i - \mathbf{r}_j|^2 } \right| \Psi \right\rangle$

${\bar f}_{\sigma \bar{\sigma}} (\gamma;\mathbf{r}) = \left \langle \Psi \left| \sum_{i=1}^{N} \delta(\mathbf{r}-\mathbf{r}_i) \sum_{j ; \sigma_i \neq \sigma_j}^{N} e^{ -\gamma |\mathbf{r}_i - \mathbf{r}_j|^2 } \right| \Psi \right \rangle$

When the wave function $\Psi$ has a standard form (sum of determinants built from molecular integrals $\phi$'s) such integrals can be easily obtained in terms of the following elementary contributions

$\phi_i(\mathbf{r}) \phi_k(\mathbf{r}) \int {\rm d} \mathbf{r}^\prime \phi_j(\mathbf{r}^\prime) \phi_l(\mathbf{r}^\prime) e^{ -\gamma |\mathbf{r} - \mathbf{r}^\prime|^2 }$

which in turn can be evaluated as generalized overlap integrals. Let us now discuss our basic approximation consisting in going from Eq. (1) to Eq. (2). This approximation can be written in a more compact way as

$\frac{ \langle \sqrt{-\log_2 f} \rangle }{ \sqrt{ -\log_2 \langle f \rangle} } \sim_{\{\gamma {\rm ~large}\}} 1 \qquad (3)$

where the symbol $\langle Q \rangle$ denotes the integration of $Q \Psi^2$ over all-particle coordinates except the $i$-th one. For a given electronic configuration $(\mathbf{r}_1,...,\mathbf{r}_N)$ and $\gamma$ large enough, $f$ is dominated by a single exponential, namely $e^{-\gamma |\mathbf{r}_i-\mathbf{r}_{j_{\rm min}}|^2}$, where $|\mathbf{r}_i-\mathbf{r}_{j_{\rm min}}|$ is the distance between the reference electron $i$ located at $\mathbf{r}$ and the closest electron labelled $j_{\rm min}$. The validity of our basic approximation is directly related to the amount of fluctuations of the quantity $f$ when various electronic configurations are considered. Note that for a given electron $j$ the distance $|\mathbf{r}_i-\mathbf{r}_j|$ can vary a lot but it is much less the case for $|\mathbf{r}_i-\mathbf{r}_{j_{\rm min}}|$ where the electron number $j_{\rm min}$ can be different from one configuration to another. When these fluctuations are small, the ratio in Eq. (3) is close to one and the approximation is of good quality. To see what happens for larger fluctuations let us write

$f = f_{\rm min} + \delta f.$

$\frac{ \langle \sqrt{-\log_2 f} \rangle }{ \sqrt{ -\log_2 \langle f \rangle} } = 1 + O[{(\delta f)}^2]$

showing that at first order in the fluctuations the ratio is still equal to one, illustrating the validity of our approximation.

A last point to discuss is the value of $\gamma$ to be chosen in practice. Because of our approximation, the limit $\gamma \rightarrow +\infty$ cannot be taken since the ratio in Eq. (3) goes to zero.[#Note] Therefore, the value of $\gamma$ has to be large enough to discriminate between the closest electron located at $r_{j_{\rm min}}$ from the other ones located at larger distances of electron $i$, while staying in the regime where the ratio in Eq. (3) stays close to one. We have found that a value of $\gamma$ depending on $\mathbf{r}$ and chosen on physical grounds allows to recover systematically the essential features of the EPLF images calculated with QMC, that is to say, with the exact expression of the {\it Min} function. To be effective, the discrimination of the closest electron with the other ones must be properly implemented. To do that, the value of $\gamma$ is adapted to keep the leading exponential $e^{-\gamma |\mathbf{r}_i-\mathbf{r}_{j_{\min}}|^2}$ significantly larger than the sub-leading exponential $e^{-\gamma |\mathbf{r}_i-\mathbf{r}_{j_{\rm next-min}}|^2}$ associated with the second closest electron $j_{\rm next-min}$. First, we define a sphere $\Omega(\mathbf{r}_i)$ centered on $\mathbf{r}_i$ with a radius $d_\Omega(\mathbf{r}_i)$. Then, locally, we represent our system made of the electron located at $\mathbf{r}_i$ and its two closest neighbors by a model system of three independent particles. If one calculates the probability of finding all the three particles inside the sphere, one finds

$P_\Omega (\mathbf{r}_i) = \left( \frac{1}{3} \int_{\Omega(\mathbf{r}_i)} {\rm d} \mathbf{r} \rho(\mathbf{r}) \right)^3$

If the density $\rho(\mathbf{r})$ is supposed constant and equal to $\rho(\mathbf{r}_i)$, the radius $d_\Omega(\mathbf{r}_i)$ of the sphere can be set such that $P_\Omega$ is equal to a fixed value

$d_\Omega(\mathbf{r}_i) = \left( \frac{4\pi}{9} {P_\Omega}^{-1/3} \rho(\mathbf{r}_i) \right)^{-1/3}$

Then, $\gamma(\mathbf{r}_i)$ is chosen in order to set a constant ratio $\kappa$ between the width of $e^{-\gamma r_{ij}^2}$ and the radius of the sphere

$\kappa = \sqrt{2 \gamma(\mathbf{r}_i)} d_\Omega(\mathbf{r}_i)$

We obtain an expression of $\gamma(\mathbf{r}_i)$ which depends on the electron density:

$\gamma(\mathbf{r}_i) = \frac{\kappa^2}{2} \left(\frac{4\pi}{9} {P_\Omega}^{-1/3} \rho(\mathbf{r}_i) \right)^{2/3}$

In our simulations, we have found that the EPLF images obtained with QMC are properly recovered using $P_\Omega = 0.001$ and $\kappa=50$.

# Note

For a given electronic configuration $(\mathbf{r}_1,...,\mathbf{r}_N)$ $f$ goes to $\exp(-\gamma |\mathbf{r}-\mathbf{r}_{\min}|^2)$ as $\gamma$ becomes large. Accordingly, in this limit $\langle \sqrt {-\log_2 f} \rangle$ goes to a finite positive constant given by $\langle |\mathbf{r}-\mathbf{r}_{\min}|^2 \rangle$. On the other hand $\langle f \rangle$ can be written as the integral of $\exp(-\gamma |\mathbf{r}-\mathbf{r}_{\min}|^2)$ over all possible configurations which can be rewritten as $\frac{1} {\sqrt{\gamma}} \sqrt{\gamma} \exp(-\gamma |\mathbf{r}-\mathbf{r}_{\min}|^2)$, that is, at large $\gamma$, the product of $1/\sqrt{\gamma}$ times a Dirac distribution. Finally, in the large-$\gamma$ limit the ratio $\frac{\langle \sqrt{-\log_2 f} \rangle }{ \sqrt{ -\log_2 \langle f \rangle} }$ goes to zero as $\sim 1/\sqrt{\gamma}$.