Extreme statistics and index distribution in the classical 1d Coulomb gas

In the last two decades, extreme value statistics in correlated random variables has received a resurgence of interest [1, 2] with the discovery of the Tracy–Widom (TW) distribution in the context of random matrix theory (RMT) [3, 4]. Since then the TW distribution has appeared ubiquitously in physics [5, 6], mathematics [7, 8] and information theory [9]. In physics it has appeared in stochastic growth models belonging to the Kardar–Parisi–Zhang (KPZ) universality class [10–16], nonintersecting Brownian motions [17], noninteracting fermions in a 1D trapping potential [18–20], disordered mesoscopic systems [21] and even in the Yang–Mills gauge theory in two dimensions [17]. It has also been measured experimentally in several systems including liquid crystals [22], coupled fiber lasers [23], or disordered superconductors [24].

Originally, the TW distribution was discovered as the limiting distribution of the largest eigenvalue $x_{\max}$ of an $N \times N$ Gaussian random matrix for which the joint probability density function (PDF) $\mathcal{P}(\{x_i\})$ of the N real eigenvalues $\{x_1, x_2, ..., x_N \}$ is known explicitly [25, 26]:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \mathcal{P}(\{x_i\}) = B_N \, {\rm e}^{-\frac{1}{2\sigma^2} \sum\nolimits_{i=1}^{N} x_i^2} \prod_{i<j} |x_i - x_j|^\beta \label{JPDF-rm} \nonumber \end{align} \tag{ 1 }$$

where B_N is the normalisation constant and $\beta=1, 2, 4$ is the Dyson index corresponding respectively to the Gaussian orthogonal, unitary and symplectic ensembles (GOE, GUE, and GSE respectively) [25, 27]. This distribution of N eigenvalues can, equivalently, be interpreted as the equilibrium Gibbs distribution, $\mathcal{P}(\{x_i\})=B_N~{\rm e}^{-\beta E(\{x_i\})}$, of a gas of charged particles with positions x_i's on a line with the energy E({x_i}) given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle E= \frac{1}{2\sigma^2 \beta} \sum_{i=1}^N x_i^2 - \frac{1}{2}\sum_{i \neq j}\log (|x_i-x_j|) . \label{H-log-gas-unscaled} \nonumber \end{align} \tag{ 2 }$$

The first term in the energy can be interpreted as the potential energy due to a confining harmonic potential, while the second term represents a logarithmic repulsion between any pair of charges. These two opposite energies compete with each other. The first term scales for large N as $\sim N\, x_{\rm typ}^2/(\sigma^2 \beta)$ where x_typ is the typical scale of the position of charges. The second term scales as N² since there are $N(N-1)$ pair of charges. Balancing the two energies lead to the fact that $x_{\rm typ} \sim \sqrt{N}$ for large N. This suggests a rescaling of the positions of the charges as $x_i \to \sigma \sqrt{\beta N}\, x_i$. In these rescaled variables, the energy is then given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle E= \frac{1}{2} \left[N\,\sum_{i=1}^N x_i^2 - \sum_{i \neq i} \log |x_i - x_j| \right], \label{H-log-gas} \nonumber \end{align} \tag{ 3 }$$

up to an unimportant constant. This system is often known as the log-gas [26]. In the large N limit, the average density $\rho_N(x)$ of these charges or the eigenvalues converges to an N-independent limiting density given by the Wigner semi-circular form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \rho_N(x)|_{N \to \infty}=\rho_\infty(x) = \frac{\sqrt{2-x^2}}{\pi}, \label{wig-sc} \nonumber \end{align} \tag{ 4 }$$

which has a finite support $x \in [-\sqrt{2}, ~\sqrt{2}]$. It turns out that the behaviour of the eigenvalues close to the soft edges $\pm \sqrt{2}$ have universal features—most easily demonstrated by the largest eigenvalue $x_{\max} = \max\nolimits_{1 \leqslant i \leqslant N} x_i$. In the log-gas picture, $x_{\max}$ corresponds to the position of the rightmost charge. Its average value is $\langle x_{\max} \rangle \sim \sqrt{2}$ for large N and it corresponds to the right edge of the semi-circle. However, $x_{\max}$ typically fluctuates from sample to sample on a scale of width N^−2/3 around the mean. The probability distribution of these typical fluctuations is described by the celebrated TW distribution. Indeed, the cumulative distribution $Q(w, N) = {\rm Prob}(x_{\max}\leqslant w, N)$, takes the scaling form for $w - \sqrt{2} = O(N^{-2/3})$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Q(w,N) \approx {{\mathcal F}}_\beta \left(\sqrt{2}N^{2/3}(w-\sqrt{2})\right) \label{TW} , \nonumber \end{align} \tag{ 5 }$$

where ${{\mathcal F}}_\beta(x)$ is the TW distribution, computed by Tracy and Widom for $\beta = 1, 2$ and 4 in terms of the solution of a Painlevé II equation [3]. For example, for $\beta = 2$ (the GUE case)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \mathcal{F}_2(x)&= \exp \left(- \int_x^\infty (y-x)q^2(y)~{\rm d}y \right),~~{\rm where},\nonumber \\ \frac{{\rm d}^2q(y)}{{\rm d}y^2}&=2 q(y)^3 + y~q(y),~{\rm with},~q(y \to \infty) \to {\rm Ai}(y), \label{TW-2} \nonumber \end{align} \tag{ 6 }$$

and ${\rm Ai}(y)$ is the Airy function. For general β, the PDF ${{\mathcal F}}'_\beta(x)$ of the TW-scaling function has non-Gaussian tails

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \mathcal{F}'_\beta(x) \sim \left\{\begin{array}{@{}lll@{}} \exp\left[-\frac{\beta}{24}|x|^3 \right]&{\rm as}&x \to - \infty \nonumber \\ \exp\left[-\frac{2\beta}{3}x^{3/2} \right]&{\rm as}&x \to \infty. \end{array}\right. \label{TW-asymptotic} \nonumber \end{align} \tag{ 7 }$$

While the typical fluctuations of $x_{\max}$ around its mean are described by the TW distribution, the atypical large fluctuations of $x_{\max}$, far from its mean to the left and right, are not described by TW but rather by the left and right large deviation tails

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{regimes_gaussian_CDF} \partial_w\,Q(w,N)\approx \left\{\begin{array}{@{}lc@{}} \exp\left[-\beta N^2 \Phi_{-}(w)\right],\,&w<\sqrt{2} \; \& \; |w-\sqrt{2}|\sim\mathcal{O}(1) \nonumber \\ \sqrt{2} N^{\frac{2}{3}}{{\mathcal F}}'_\beta \left(\sqrt{2} N^{\frac{2}{3}}(w-\sqrt{2})\right), & |w-\sqrt{2}|\sim\mathcal{O}(N^{-\frac{2}{3}}) \nonumber \\ \exp\left[-\beta N \Phi_{+}(w)\right], & w>\sqrt{2} \, \& \, |w-\sqrt{2}|\sim\mathcal{O}(1) , \end{array}\right. \nonumber \end{align} \tag{ 8 }$$

where the right large deviation function (LDF) $\Phi_+(w)$ was obtained explicitly for $\beta = 1$ in [28] and for arbitrary β in [31]. On the other hand, $\Phi_-(w)$ was computed explicitly for all β in [29, 30]. It was argued that the left and the right large deviation tails can be interpreted as the free energies of two different thermodynamic phases of the Coulomb gas, separated by a third order phase transition [33] in the large N limit. Similar third order phase transitions have also been found in a variety of other systems [32–37], including in higher dimensions $\newcommand{\g}{g} d\geqslant 1$ [38–41].

The TW distribution was initially derived for an harmonic potential. However, it was found later that the typical distribution of $x_{\max}$ is universally given by the TW distribution, irrespective of the shape of the confining potential. This holds provided the average charge/eigenvalue density has a finite support and moreover, vanishes as a square root at the upper edge of the support (for a recent review see [48]). One then naturally asks the question: what happens to the universality of the TW distribution, if instead of the confining potential, one changes the form of the repulsive pairwise interaction? A natural setting to address this question corresponds to a model in 1d of N charged particles in presence of a confining harmonic potential and just replacing the logarithmic pairwise repulsion by the true Coulomb repulsion in 1d, i.e. a linear $|x_i - x_j|$ interaction term in equation (3) instead of $\log |x_i - x_j|$.

It turns out that the resulting model is not purely academic. Indeed, this is a model that has been well studied in the physics literature in the context of 1d charged plasma [49]—known as the one-dimensional one component plasma (1d OCP) or the 'jellium' model. In the jellium model, there are two species of opposite charges with a vanishing total charge. Each pair of charges interact with each other via the Coulomb potential which happens to be linear in 1d. If one treats one of the species, say the ones with negative charges, in a mean-field fashion with a uniform background density, this generates an effective confining harmonic potential for the positive charges. In addition, the positive charges repel each other via the linear Coulomb interaction. This gives rise precisely to an energy as in equation (3) with the logarithm replaced by a linear term. This model is a paradigm for 1d charged plasma [49] as several observables can be calculated analytically [50–54]. For this model most of the earlier studies considered bulk properties in the thermodynamic limit. In a recent letter [42], we addressed the extreme value question in the 1d OCP or the 'jellium' model where we showed analytically that the limiting distribution of the typical fluctuations of $x_{\max}$ is indeed different from the TW distribution. Moreover, by computing the left and the right LDFs explicitly, we have shown that the third-order phase transition between a pushed gas (left large deviation) and a pulled gas (right large deviation) is still present in this system as in the case of the log-gas. One of the purposes of the current paper is to provide a detailed derivation of these results presented in the letter [42].

In fact, the question of universality with respect to the pairwise repulsion term is not restricted just to the rightmost particle position $x_{\max}$, but can also be addressed for other observables. For instance, one can ask how sensitive is the statistics of the gap g between the positions of the rightmost and the next rightmost particles for large N, as one changes the form of the pairwise interaction? Indeed, for the GUE, the PDF of the gap takes the scaling form, for large N

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_G(g,N) \approx \sqrt{2}N^{2/3}~h_2(g \sqrt{2}N^{2/3}) , \label{gap_GUE} \nonumber \end{align} \tag{ 9 }$$

where the scaling function h₂(x) was computed explicitly in [43, 44]. In this paper, we compute exactly the gap distribution in the 'jellium' model and show that it is different from the GUE-log-gas in equation (9).

Another interesting observable that has been studied extensively in the context of random matrices is the index N₊  that denotes the number of positive eigenvalues, or equivalently the number of charges on the positive semi-axis. Obviously N₊ is a random variable with values $0 \leqslant N_+ \leqslant N$. It was shown that N₊ typically fluctuates around its mean value N/2 on a scale of width $\sqrt{\log N}$ and the typical fluctuations are given by a Gaussian form [45, 46, 47]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle P_I(N_+,N) \approx \exp \left[- \frac{\beta \pi^2}{2 \ln N}(N_+-N/2)^2\right]. \label{gauss_index} \nonumber \end{align} \tag{ 10 }$$

The atypical large deviations of (N₊  −  N/2)  =  O(N) were computed for large N and it was found that [46, 47]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{index_ldf} P_I(N_+,N) \sim {\rm e}^{- \beta N^2 \Psi\left(\frac{N_+}{N}\right)} \nonumber \end{align} \tag{ 11 }$$

where the rate function $\Psi(c)$ has a logarithmic singularity at c  =  1/2. In this paper, we study this index distribution analytically for the 'jellium' model and find that it is rather different from the log-gas case.

Thus the purpose of this paper is essentially twofold:

• to present the detailed calculations of the distribution of $x_{\max}$ in the 'jellium' model
• to present new exact results for the distribution of two other observables in the 'jellium' model: (i) the gap g between the positions of the two rightmost particles and (ii) the index N₊ denoting the number of particles on the positive semi-axis.

Interestingly, as we will show, the function that characterises the limiting distribution of $x_{\max}$ happens to also characterise the limiting distribution of the gap and that of the index. One of the main results of this paper is to show that changing the pairwise interaction indeed changes significantly the behaviour of these observables, thereby changing their universality class.

We consider N charges on a line, with positions {x_i}, confined by an external harmonic potential and interacting pairwise via the true 1d repulsive Coulomb potential. We assume that the system is in thermal equilibrium such that the probability to observe the system in a configuration {x_i} is given by the Boltzmann distribution

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \mathcal{P}(x_1,x_2,...,x_N) = \frac{\exp \left[-\beta\,E(x_1,x_2,...,x_N)\right]}{Z_N}, \label{jpdf-1dc} \nonumber \end{align} \tag{ 12 }$$

where β is the inverse temperature, Z_N is the normalization constant and the energy of the configuration is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle E(\{x_i\})=\frac{N}{2} \sum_{i=1}^N x_i^2 - \alpha \sum_{i \neq j} |x_i-x_j| , \label{E-1dc} \nonumber \end{align} \tag{ 13 }$$

where $\newcommand{\g}{g} \newcommand{\al}{\alpha} \alpha \geqslant 0$ denotes the strength of the Coulomb repulsion. As in the log-gas case in equation (3), the prefactor N in the first term ensures that the x_i's are of order $O(1)$.

Let us first consider two limiting temperature regimes: (i) very high temperature when $\beta \ll N$ and (ii) very low temperature with $\newcommand{\g}{g} \beta \gg N$. In the former case, the interaction between the charges become totally irrelevant and the particles behave as N independent random variables with Gaussian distributions. In contrast, in case (ii) the interaction term dominates and the positions of the particles get 'frozen' at equidistant points in the interval $\newcommand{\al}{\alpha} [-2\alpha, + 2 \alpha]$, with very small fluctuations around them. While one can study this model at all temperatures exactly (see for example the results for $x_{\max}$ at very high and very low temperature in appendix), it turns out that the most interesting situation occurs when $\beta = O(N)$ when both the interaction term as well as the confining potential compete with each other. Henceforth, in this paper, we will focus on the case where $\beta = N$ so that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \beta E(\{x_i\})=\frac{N^2}{2} \sum_{i=1}^N x_i^2 - \alpha N\, \sum_{i \neq j} |x_i-x_j| . \label{E-1dc_2} \nonumber \end{align} \tag{ 14 }$$

This is the so-called 'jellium' model, whose bulk properties in the thermodynamic limit ($N \to \infty$) have been studied extensively before [50–54]. In this paper, our focus is on the edge behaviour, in particular the distribution of the position of the rightmost particle $x_{\max}$, as well as the gap g between the positions of the two rightmost charges. In addition, we also compute the distribution of the index N₊ in the large N limit. Let us summarise our main results:

• (a)  
  Distribution of $x_{\max}$: it is well known that, in the large N limit, the average density of charges (normalised to unity) is uniform $\newcommand{\f}{f} \newcommand{\al}{\alpha} \rho_N(x)|_{N\to \infty} = \rho_\infty(x) = \frac{1}{4\alpha}$, for $\newcommand{\al}{\alpha} -2\alpha \leqslant x \leqslant 2 \alpha$. Thus the average $\newcommand{\al}{\alpha} \langle x_{\max}\rangle \to 2 \alpha$. For large but finite N, $x_{\max}$ fluctuates around this mean value with typical fluctuations scaling as $O(1/N)$. Indeed we compute the full cumulative distribution

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Q(w,N)={\rm Prob}.\,[x_{\max}\leqslant w, N], \label{Q-1dc} \nonumber \end{align} \tag{ 15 }$$

  and show that it exhibits three different regimes (see figure 1)

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \displaystyle Q(w,N) \approx \left\{\begin{array}{@{}ll@{}} \ee^{- N^3\, \Phi_{-}(w) + O(N^2)},&0< 2\al-w \sim O(1) \nonumber \\ F_\al(N(w-2\alpha) + 2 \alpha),&|w-2\al| \sim O(1/N) \nonumber \\ 1-\ee^{- N^2\, \Phi_+(w)+O(N)},& 0<w-2\al \sim O(1) . \end{array}\right. \label{result-1} \nonumber \end{align} \tag{ 16 }$$

  The second line denotes the regime for typical fluctuations $\newcommand{\al}{\alpha} |w - 2 \alpha| \sim O(1/N)$ where the scaling function $\newcommand{\al}{\alpha} F_\alpha(x)$ satisfies a nonlocal eigenvalue equation

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle \frac{{\rm d}F_\alpha(x)}{{\rm d}x}= A(\al)\, \ee^{-x^2/2}\, F_\alpha(x+4\al) , \label{limit.00} \nonumber \end{align} \tag{ 17 }$$

  with $\newcommand{\al}{\alpha} A(\al)$ as the unique eigenvalue that can be determined (see later). The tails of the distribution $\newcommand{\al}{\alpha} F_\al(x)$ are given by

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle F_\al'(x) \approx \left\{\begin{array}{@{}lll@{}} \exp\left[-|x|^3/{24\al}+ O(x^2)\right] & {\rm as} & x \to -\infty \nonumber \\ \exp\left[-x^2/2 + O(x)\right] & {\rm as} & x \to \infty . \end{array}\right. \label{fx_asymp-0} \nonumber \end{align} \tag{ 18 }$$

  This is the analogue of the TW distribution found in the log-gas case (see equation (8)). Clearly, this limiting distribution is different from the TW distribution.The first and third lines in equation (16) describe the atypical large fluctuations of $x_{\max}$, respectively to the left and the right of the central typical regime. We compute explicitly the rate functions $\Phi_-(x)$ and $\Phi_+(x)$:

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{sm-left_rf-0} \Phi_{-}(w) = \left\{\begin{array}{@{}lc@{}} \frac{(2\al-w)^3}{24\al}& -2\al\leqslant w \leqslant 2\al \nonumber \\ \frac{w^2}{2} + \frac{2}{3} \al^2, & w \leqslant -2\al. \end{array}\right. \nonumber \end{align} \tag{ 19 }$$

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{sm-phi+0} \Phi_+(w) = \frac{(w-2\al)^2}{2} , \quad \;\;\;\quad w > 2\al . \nonumber \end{align} \tag{ 20 }$$

  The three regimes in equation (16) are shown schematically in figure 1. It is easy to check that the central part described by $\newcommand{\al}{\alpha} F_\al(x)$ in equation (16) matches smoothly with the two large deviation regimes flanking this central part.
• (b)  
  Distribution of the gap: we compute the distribution P_G(g, N) of the gap g between the positions of the two rightmost particles. We show that, for large N, it has the scaling form

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle \label{eq:scaling_g} P_G(g,N) \approx N~h_\alpha(gN) , \nonumber \end{align} \tag{ 21 }$$

  where the scaling function $\newcommand{\al}{\alpha} h_\alpha(z)$ is given by

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle \label{eq:halpha} h_\alpha(z) = \Theta(z)\,A(\alpha) \int_{-\infty}^\infty {\rm d}y~(y+z-4\alpha)~{\rm e}^{-(y+z-4\alpha)^2/2}~F_\alpha(y) . \nonumber \end{align} \tag{ 22 }$$

  In equation (22), $\newcommand{\al}{\alpha} F_\alpha(x)$ is again the unique solution of equation (17) with eigenvalue $\newcommand{\al}{\alpha} A(\alpha)$ and $\Theta(z)$ is the Heaviside step function.
• (c)  
  Distribution of the index: we have also computed analytically, for large N, the distribution $P_I(N_+, N)$ of the index N₊ , i.e. the number of charges $N_+= \sum\nolimits_{i=1}^N \theta(x_i)$ on the positive semi-axis. From the symmetry of the energy E({x_i}) about the origin in (13), it is evident that $\langle N_+ \rangle =N/2$ and furthermore the full distribution $P_I(N_+, N)$ is symmetric around N₊  =  N/2. Indeed, for large N, we show that it approaches a scaling form

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle P_I(N_+, N) \approx 4\alpha\, f_\alpha\left(4 \alpha \left(N_+-{N}/{2}\right)\right)\, , \label{limiting.1_intro} \nonumber \end{align} \tag{ 23 }$$

  where the scaling function $\newcommand{\al}{\alpha} f_\alpha(z)$ is given by

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle f_\alpha(z)= \frac{F_\alpha(z+2\alpha)\, F_\alpha(-z+2\alpha)}{\int_{-\infty}^{\infty} {\rm d}z\, F_\alpha(z+2\alpha)\, F_\alpha(-z+2\alpha)} \; \label{limiting.2_intro} \nonumber \end{align} \tag{ 24 }$$

  where $\newcommand{\al}{\alpha} F_\alpha(x)$ is again the unique solution of (17). Note that the typical fluctuations of N₊ around its mean are here of order $O(1)$ (see equation (23)), while they are of order $O(\sqrt{\ln N})$ in the log-gas (see equation (10)). Moreover, we show that the scaling function $\newcommand{\al}{\alpha} f_\alpha(z)$ has non-Gaussian tails

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle f_\alpha(z) \sim \exp\left[- \frac{1}{24\alpha}\,|z|^3\right]\quad {\rm as}\quad |z|\to \infty \, . \label{asymp.1_intro} \nonumber \end{align} \tag{ 25 }$$

  This limiting distribution is thus non-Gaussian, unlike in the log-gas case where it is known to be Gaussian (see equation (10)). The function $\newcommand{\al}{\alpha} f_\alpha(z)$ describes only the typical fluctuations of N₊ of order $O(1)$ around its mean. The atypical fluctuations of order $O(N)$ on both sides of the mean are described by symmetric large deviation tails:

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle P_I(N_+=cN,N)\approx \exp \left(-N^3 \, \Psi(c) \right) \;, \;\; \Psi(c) = \frac{8 \alpha^2}{3}~|c-1/2|^3 ,~~0\leqslant c \leqslant 1 . \label{LDF1} \nonumber \end{align} \tag{ 26 }$$

  The rate function $\Psi(c)$ is simple here and is rather different from the corresponding one in the log-gas case [46, 47].

Zoom In Zoom Out Reset image size

With the choice $\beta = N$, the partition function of the model is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle Z_N = \int_{-\infty}^\infty {\rm e}^{- \beta E(\left\{x_i\right\})} {\rm d}x_1 \cdots {\rm d}x_N , \;{\rm with}\;\quad \beta E(\left\{x_i\right\}) = \frac{N^2}{2} \sum_{i=1}^N x_i^2 - \alpha N\, \sum_{i \neq j} |x_i-x_j| . \label{partition_func} \nonumber \end{align} \tag{ 27 }$$

In the large N limit (equivalently the zero temperature limit), the partition function is dominated by the ground state (minimum energy configuration). To find this minimum energy configuration it is convenient to first rewrite the partition function Z_N in (27) using the fact that E({x_i}) is symmetric under permutations of the x_i's. Hence

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z_N = N ! \int_{-\infty}^\infty \cdots \int_{-\infty}^\infty {\rm d}x_1 \cdots {\rm d}x_N \; {\rm e}^{- \beta E(\left\{x_i\right\})} \prod_{j=2}^N \Theta(x_j - x_{j-1}) . \label{partition_func2} \nonumber \end{align} \tag{ 28 }$$

For such an ordered configuration $x_1<x_2<\cdots<x_N$, we can eliminate the absolute values and rewrite the energy function as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \beta E(\{x_i\})&=\frac{N^2}{2} \sum_{i=1}^N x_i^2 - 2\alpha N \sum_{i > j} (x_i-x_j), \nonumber \\ &=\frac{N^2}{2} \sum_{i=1}^N x_i^2 - 2\alpha N \sum_{i=1}^N (2i-N-1)x_i, \nonumber \\ &=\frac{N^2}{2} \sum_{i=1}^N \left(x_i^2 - \frac{4\alpha}{N}(2i-N-1)x_i \right), \nonumber \\ &= \frac{N^2}{2} \sum_{i=1}^N \left(x_i - \frac{2 \alpha}{N}(2i-N-1)\right)^2 + C_N(\alpha), \label{E-ordered} \nonumber \end{align} \tag{ 29 }$$

where $\newcommand{\al}{\alpha} C_N(\alpha)= 2 \alpha^2 \sum\nolimits_{i=1}^N \left(2i-N-1\right){}^2$ is just a constant. The minimum energy configuration corresponds to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle x_i = x_i^*=\frac{2 \alpha}{N}(2i-N-1),~~{\rm for}~i=1,2,...,N.\label{x_j-min} \nonumber \end{align} \tag{ 30 }$$

This implies that, in the minimum energy configuration the charges are placed at regular intervals of length $\newcommand{\f}{f} \newcommand{\al}{\alpha} \frac{4 \alpha}{N}$. The rightmost particle is at $\newcommand{\al}{\alpha} x_N^* = 2 \alpha(1-1/N)$ and the leftmost particle is at the symmetrically opposite place $\newcommand{\al}{\alpha} x_1^* = -2 \alpha(1-1/N)$. Hence it is clear that the charge density is supported over a finite support and in the large N limit it is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \rho_\infty(x) = \left\{\begin{array}{@{}ccc@{}} \frac{1}{4\alpha},&{\rm for}&-2\alpha \leqslant x \leqslant 2 \alpha \nonumber \\ 0& &{\rm otherwise} \end{array}\right. . \label{rho-1dc} \nonumber \end{align} \tag{ 31 }$$

Note that $\rho_\infty(x)$ is different from the Wigner semi-circle (4) obtained in the log-gas.

For this jellium model, different thermodynamic properties have been studied extensively [49–54]. In particular, Baxter [52] analysed the partition function Z_N,L of the jellium model confined in a finite box $[-L, L]$, i.e. the following multiple integral

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:ZNL} Z_{N,L} = N ! \int_{-L}^L \cdots \int_{-L}^L {\rm d}x_1 \cdots {\rm d}x_N \; {\rm e}^{- \beta E(\left\{x_i\right\})} \prod_{j=2}^N \Theta(x_j - x_{j-1}) . \nonumber \end{align} \tag{ 32 }$$

In computing this integral (32), Baxter introduced [52], as an intermediate step, an auxiliary function $\newcommand{\al}{\alpha} F_\alpha(x)$ that satisfies a non-local eigenvalue equation defined in equation (17). However, the observables in this model that we are interested in have not been studied, to the best of our knowledge. Remarkably, we find that the same auxiliary $\newcommand{\al}{\alpha} F_\alpha(x)$ function that Baxter introduced for the analysis of the partition function of the system in a finite box, also plays a key role in determining the distributions of our observables on the infinite line.

In this section, we focus on the position $x_{\max}$ of the rightmost particle on the infinite line. From the analysis of the average density in equations (30) and (31), it is clear that, in the limit $N \to \infty$, the mean position of the rightmost particle is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle \label{av_xmax} \langle x_{\max}\rangle = x_N^* \approx 2 \alpha . \nonumber \end{align} \tag{ 33 }$$

To derive the distribution of $x_{\max}$, it is convenient to consider the cumulative distribution

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle Q(w,N) = {\rm Prob}.[x_{\max} \leqslant w] = {\rm Prob.}(x_1\leqslant w, \cdots,x_N\leqslant w).\nonumber \end{align*}$$

Using the Boltzmann distribution $\mathcal{P}(x_1, x_2, ..., x_N)$ in equation (12), one can express $Q(w, N)$ as the ratio of two partition functions

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle Q(w,N)= \frac{Z_N(w)}{Z_N(\infty)},~{\rm where,}\label{Q_N-1} \nonumber \end{align} \tag{ 34 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ee}{{\rm e}} \displaystyle Z_N(w)= \int_{-\infty}^w {\rm d}x_1\cdots \int_{-\infty}^w {\rm d}x_N~\ee^{-\beta E(\{x_i\})}, \label{Z_N-w} \nonumber \end{align} \tag{ 35 }$$

where $\beta\, E(\{x_i\})$ is given in (13) and $\newcommand{\e}{{\rm e}} Z_N(\infty) \equiv Z_N$ given in equation (27). Again, it is convenient to work with ordered configurations of the x_i's, $-\infty < x_1\leqslant x_2 \leqslant...\leqslant x_N \leqslant w$ as before and one gets (using equation (29))

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle Z_N(w) \propto \int_{-\infty}^w {\rm d}x_N \int_{-\infty}^{x_N} {\rm d}x_{N-1}...\int_{-\infty}^{x_2} {\rm d} x_{1}~\ee^{- \frac{N^2}{2} \sum\nolimits_{i=1}^N \left(x_i - \frac{2 \alpha}{N}(2i-N-1)\right)^2} \label{Z_N-w-ord} \nonumber \end{align} \tag{ 36 }$$

where we have replaced the product of theta functions by constraining the limits of the integrals. It is natural now to make a change of variables

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \displaystyle \ep_i=\left[Nx_i-2\al(2i-N-1)\right],~i=1,2,\cdots,N . \label{transformation} \nonumber \end{align} \tag{ 37 }$$

The ordering condition $x_{i-1} < x_i$ translates to the following constraint on the $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \epsilon_i$'s

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \displaystyle \label{cond_epsilon} \ep_{i-1}<\ep_i+4 \alpha ,~{\rm for}~i=2,3,...,N . \nonumber \end{align} \tag{ 38 }$$

The last constraint x_N  <  w in equation (36) translates to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \displaystyle \label{cond_epsilon_N} \epsilon_N < N\,(w-2\al) + 2 \alpha . \nonumber \end{align} \tag{ 39 }$$

Consequently, Z_N(w) reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle Z_N(w) \propto D_\alpha\left(N(w-2\al)+2\al,N\right), \label{Z_N_D_a} \nonumber \end{align} \tag{ 40 }$$

where the function $\newcommand{\al}{\alpha} D_\al(x, N)$ on the right hand side (rhs) is given by the N-fold integral

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle D_\al(x,N) =~\int_{-\infty}^{x} {\rm d} \ep_N\int_{-\infty}^{\ep_N+4 \al}{\rm d} \ep_{N-1}\dots \int_{-\infty}^{\ep_2+4\al} {\rm d}\ep_1~\ee^{-\frac{1}{2}\sum\nolimits_{i=1}^N \ep_i^2} . \label{D_N} \nonumber \end{align} \tag{ 41 }$$

We remark that in the original jellium model in equation (27), the interaction between the x_i's is long-ranged (as every charge is coupled to every other charge). Remarkably however, after the ordering of the positions and the change of variables in equation (37), the interactions between the new variables $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \epsilon_i$'s in equation (41) become short-ranged, i.e. $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \epsilon_i$ interacts only with its two nearest neighbours $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \epsilon_{i-1}$ and $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \epsilon_{i+1}$. Note that this however is true only for $\newcommand{\al}{\alpha} \alpha > 0$ strictly (for $\newcommand{\al}{\alpha} \alpha = 0$ the gas is non-interacting). Therefore, the function $\newcommand{\al}{\alpha} D_\al(x, N)$ in equation (41) can be interpreted as a restricted partition function of this constrained short-ranged interacting gas.

To make further progress, we substitute Z_N(w) from equations (40) and (41) into equation (34) and obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle Q(w,N) = \frac{D_\al(x,N)}{D_\al(\infty,N)} \equiv F_\alpha(x,N) . \label{def_FN} \nonumber \end{align} \tag{ 42 }$$

Taking derivative with respect to x in (42), and using (41), we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{dF_N.1} \frac{{\rm d}F_\al(x,N)}{{\rm d}x} = \frac{D_{\al}(\infty,N-1)}{D_{\al}(\infty,N)}\,\ee^{-\frac{x^2}{2}} F_{\al}(x+4\al,N-1) . \nonumber \end{align} \tag{ 43 }$$

These equations (42) and (43) are actually exact for all x and N. To make progress, we will consider the $N \to \infty$ limit. In this limit, the typical fluctuations of $x_{\max}$ around its mean value $\newcommand{\al}{\alpha} 2 \alpha$ turn out to be of order $O(1/N)$, while atypical large fluctuations can be of order $O(1)$. Below, we analyse the probability distribution of typical and atypical fluctuations separately.

4.1. Typical fluctuations of ${x_{\max}}$

To analyse the typical fluctuations, we need to keep the argument $\newcommand{\al}{\alpha} x = (w-2\al) N + 2 \al$ of $\newcommand{\al}{\alpha} D_\alpha(x, N)$ fixed in equation (40), while we take the $N \to \infty$ limit. In addition, we need to estimate the ratio $\newcommand{\f}{f} \newcommand{\al}{\alpha} \frac{D_{\al}(\infty, N-1)}{D_{\al}(\infty, N)}$ in equation (43) in the large N limit. As discussed earlier, since $\newcommand{\al}{\alpha} D_\al(\infty, N)$ is the partition function of a short-ranged gas, we expect that its free energy $\newcommand{\al}{\alpha} -\ln D_\al(\infty,N)$ is extensive in N. Hence it follows that, for $\newcommand{\al}{\alpha} \alpha >0$, $\newcommand{\al}{\alpha} D_\al(\infty, N) \sim [A(\alpha)]^{-N}$ for large N, where $\newcommand{\al}{\alpha} \ln A(\alpha)$ is the free energy per particle of the short-ranged interacting gas. Hence, for $\newcommand{\al}{\alpha} \alpha >0$, the ratio

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \frac{D_{\al}(\infty,N-1)}{D_{\al}(\infty,N)} \to A(\alpha),~{\rm as}~N \to \infty . \label{D-inf} \nonumber \end{align} \tag{ 44 }$$

On the other hand, $\newcommand{\al}{\alpha} \alpha = 0$ is different. In this case, evaluating the integral in equation (41) exactly leads to $D_0(\infty, N) =(2\pi){}^{N/2}/N!$. Hence the ratio

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{ratio_a0} \frac{D_0(\infty,N-1)}{D_0(\infty,N)} = \frac{N}{\sqrt{2 \pi}} . \nonumber \end{align} \tag{ 45 }$$

In the discussion below, we will restrict ourselves to the $\newcommand{\al}{\alpha} \alpha > 0$ case. Substituting the result in equation (44) on the rhs of equation (43) and anticipating further that the function $\newcommand{\al}{\alpha} F_\al(x, N)$ converges to a limiting form $\newcommand{\al}{\alpha} F_\al(x)$ for large N, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle F_\alpha(x,N \to \infty) = F_\alpha(x), \label{limit.0} \nonumber \end{align} \tag{ 46 }$$

we find that $\newcommand{\al}{\alpha} F_\al(x)$ satisfies a nonlocal equation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle \frac{{\rm d}F_\alpha(x)}{{\rm d}x}= A(\al)\, \ee^{-x^2/2}\, F_\alpha(x+4\al) . \label{limit.1} \nonumber \end{align} \tag{ 47 }$$

The prefactor $\newcommand{\al}{\alpha} A(\al)$ on the rhs is still unknown. The function $\newcommand{\al}{\alpha} F_\alpha(x)$ is a cumulative probability distribution and hence satisfies the positivity condition $\newcommand{\al}{\alpha} 0\leqslant F_\alpha(x)\leqslant 1$ for $-\infty < x < \infty$, along with the boundary conditions $\newcommand{\al}{\alpha} F_\alpha(-\infty)=0$ and $\newcommand{\al}{\alpha} F_\al(\infty)=1$. It turns out that the solution of equation (47) satisfies these conditions only for a specific value $\newcommand{\al}{\alpha} A(\al)$—in this sense equation (47) can be interpreted as a non-local eigenvalue equation.

As we have remarked earlier, the same non-local eigenvalue equation (47) also appeared in Baxter's analysis of Z_N,L in equation (32). Indeed, we remark that we can give a probabilistic interpretation to the integral Z_N,L in equation (32). Up to a prefactor, this is just the probability that all the particles in an infinite system are contained in $[-L, +L]$, which in turn, is the probability that the maximum of |x_i|'s is less than L, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Z_NL_prob} Z_{N,L} \propto {\rm Prob.}\, \left[\max\{|x_1|, |x_2|, \cdots, |x_N|\} < L\right]. \nonumber \end{align} \tag{ 48 }$$

While computing $\newcommand{\al}{\alpha} A(\alpha)$ analytically for all $\newcommand{\al}{\alpha} \alpha>0$ seems hard, it is possible to determine its small and large $\newcommand{\al}{\alpha} \al$ behaviours. We first consider the large $\newcommand{\al}{\alpha} \al$ limit. In this case, using the boundary condition that $\newcommand{\al}{\alpha} F_\al(x \to \infty) = 1$ on the rhs of equation (47), one obtains, as $\newcommand{\al}{\alpha} \al \to \infty$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \frac{{\rm d} F_\alpha(x)}{{\rm d}x} \approx A(\alpha) {\rm e}^{-x^2/2} \label{largeal_1} . \nonumber \end{align} \tag{ 49 }$$

Integrating over x one finds $\newcommand{\al}{\alpha} F_\al(x) \approx A(\al)\, \int_0^x {\rm e}^{-y^2/2}\, {\rm d}y$. Using once again that $\newcommand{\al}{\alpha} F_\al(x \to \infty) = 1$ gives $\newcommand{\al}{\alpha} A(\al) \approx 1/\sqrt{2 \pi}$ as $\newcommand{\al}{\alpha} \al \to \infty$. In contrast, the $\newcommand{\al}{\alpha} \alpha \to 0$ limit is less trivial. However, this was already determined by Baxter [52]. Translating his asymptotic results to our case (his notations are quite different from ours) we finally get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle A(\al) \approx \left\{\begin{array}{@{}lcl@{}} 1/(4\,{\rm e}\,\al),&{\rm as}&\al \to 0 \nonumber \\ 1/\sqrt{2 \pi},&{\rm as}&\al \to \infty. \end{array}\right. \label{a-al-asymp} \nonumber \end{align} \tag{ 50 }$$

Note that the fact that $\newcommand{\al}{\alpha} A(\alpha)$ diverges as $\newcommand{\al}{\alpha} \alpha \to 0$ is not surprising since the ratio in equation (45) for $\newcommand{\al}{\alpha} \alpha = 0$ diverges in the $N \to \infty$ limit. For other values of α, $\newcommand{\al}{\alpha} A(\alpha)$ can be computed numerically using a shooting method, as discussed later.

Asymptotic behaviours of $\newcommand{\al}{\alpha} F_\alpha(x)$. The tails of $\newcommand{\al}{\alpha} F_\alpha(x)$ to leading order can be determined for arbitrary $\newcommand{\al}{\alpha} \alpha > 0$ as it does not require the explicit knowledge of $\newcommand{\al}{\alpha} A(\alpha)$. We first consider the $x \to \infty$ limit. In this limit we replace $\newcommand{\al}{\alpha} F_\alpha(x + 4 \alpha)$ by 1 on the rhs of (47). This gives, to leading order, the Gaussian tail $\newcommand{\ee}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\e}{{\rm e}} F_\alpha'(x) \sim \ee^{-x^2/2}$ for large x. To compute the left tail, $x \to -\infty$, we make the following ansatz

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \displaystyle \label{Falpha-inf} F_\alpha(x) \approx \ee^{-a_0\,|x|^\delta},~{\rm as}~x \to -\infty, \nonumber \end{align} \tag{ 51 }$$

where a₀ and δ are to be determined. We substitute this ansatz on both sides of (47) and then equate the powers of $|x|$ in the exponential. The rhs yields $\newcommand{\ee}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\e}{{\rm e}} {\rm rhs} \approx A(\al)\, \ee^{-x^2/2-a_0\, (|x|-4\al){}^\delta}$. For large $|x|$, $\newcommand{\al}{\alpha} (|x|-4\al){}^\delta \sim |x|^\delta(1 - 4\, \alpha\, \delta/|x|)$ to leading orders. Hence the rhs behaves as

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \displaystyle {\rm rhs} \approx A(\al) \ee^{-a_0\,|x|^\delta - x^2/2 + 4\,\al\,\delta\,|x|^{\delta-1}}. \nonumber \end{align*}$$

The left hand side (lhs) of (52) behaves as

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\ee}{{\rm e}} \displaystyle {\rm lhs}\approx \ee^{-a_0\, |x|^\delta} \nonumber \end{align*}$$

to leading order. Comparing both sides, we see that the term x²/2 and $|x|^{\delta-1}$ on the rhs must cancel each other. This implies that $\delta = 3$ and $\newcommand{\al}{\alpha} a_0 = 1/(8\alpha \, \delta) = 1/(24\, \alpha)$. This provides the leading left tail $\newcommand{\al}{\alpha} F_\alpha'(x)$ in (52). Together, the leading order asymptotic tails are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle F_\al'(x) \approx \left\{\begin{array}{@{}lll@{}} \exp\left[-|x|^3/{24\al}+ O(x^2)\right] & {\rm as}& x \to -\infty \nonumber \\ \exp\left[-x^2/2 + O(x)\right] & {\rm as} & x \to \infty . \end{array}\right. \label{fx_asymp} \nonumber \end{align} \tag{ 52 }$$

Note that the leading left tail of $\newcommand{\al}{\alpha} F'_\al(x)$ is similar to the left tail of the TW distribution in (7) (with $\newcommand{\al}{\alpha} \beta = 1/\alpha$), while the right tail in equation (52) is different from the right tail in equation (7).

For general $\newcommand{\al}{\alpha} \al>0$, it is difficult to determine the eigenvalue $\newcommand{\al}{\alpha} A(\alpha)$ as well as the full scaling function $\newcommand{\al}{\alpha} F_\alpha(x)$ explicitly. However they can be obtained by solving (47) numerically by tuning the value of $\newcommand{\al}{\alpha} A(\al)$ using the standard shooting method⁵. This gives $\newcommand{\al}{\alpha} F_\al(x)$ and $\newcommand{\al}{\alpha} A(\al)$ simultaneously. In figure 2 (left panel), we plot $\newcommand{\al}{\alpha} A(\al)$ versus $\newcommand{\al}{\alpha} \al$ and compare with its predicted asymptotics in (50). In figure 2 (right panel), we compare $\newcommand{\al}{\alpha} F_\al'(x)$ evaluated numerically using this shooting method, with the one obtained from direct Monte-Carlo simulation of the jellium model. The agreement is excellent.

Zoom In Zoom Out Reset image size

4.2. Atypical large fluctuations of ${x_{\max}}$

In the previous section we have studied the typical fluctuations of $x_{\max}$ on a scale of order $O(1/N)$ around its mean $\newcommand{\al}{\alpha} 2\al$ in the large N limit. We have shown that this centred and scaled limiting cumulative distribution is described by $\newcommand{\al}{\alpha} {\rm Prob.}(x_{\max}<w) = Q(w, N) =$$\newcommand{\al}{\alpha} \approx F_\alpha(N(w-2\al) + 2\al)$ where the scaling function $\newcommand{\al}{\alpha} F_\alpha(x)$ is given in equation (47), along with the tails given in equation (52). However this limiting distribution does not describe large fluctuations of $O(1)$ at far left or right of the mean. In the log-gas case, these large deviation functions were computed exactly [33] as described in the introduction, that revealed an interesting third order phase transition between a 'pushed' and a 'pulled' phase. It is then interesting to ask whether a similar phase transition also exists in the jellium model. This motivated us to study the probability of large deviations in the jellium model. Our exact computations show that a similar third order phase transition also exists in this case. Below, we discuss the left and right large deviation functions separately as they correspond to different physics.

4.2.1. Left large deviation.

We start with equations (34) and (35), with the energy $\beta E(\{x_i\})$ given in equation (27). We need to compute the leading behaviour of the partition function Z_N(w) for large N with a wall at w such that $\newcommand{\al}{\alpha} 0 < (2 \alpha - w) \sim 1$. This can be performed as follows: One first introduces a macroscopic empirical charge density in $(-\infty, w]$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{sm:rho_w} {\rho}_w(x)= \frac{1}{N} \sum_{i=1}^N \delta(x-x_i) . \nonumber \end{align} \tag{ 53 }$$

Note that $\rho_w(x)$ is normalised to unity. In terms of $\rho_w(x)$, the energy function $\beta E(\{x_i \})$ in equation (27) can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \beta\,E[\{x_i\}] \equiv \mathcal{E}[\rho_w(x)] \nonumber \end{align} \tag{ 54 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{sm-E-Ncube} \mathcal{E}[\rho_w(x)]=N^3\,\left(\frac{1}{2} \int_{-\infty}^w {\rm d}x~x^2\rho_w(x) - \al\int_{-\infty}^w {\rm d}x\int_{-\infty}^w {\rm d}y~\rho_w(x)\rho_w(y)~|x-y|\right) . \nonumber \end{align} \tag{ 55 }$$

The N-fold integration in the partition function Z_N(w) in (35) is carried out in two steps. In the first step, we fix the macroscopic density $\rho_w(x)$ and then sum over all the microscopic configurations of x_i's consistent with this density $\rho_w(x)$. In the second step, we sum over all possible macroscopic densities $\rho_w(x)$ that are positive and normalised to unity $\int_{-\infty}^w \rho_w(x) \, {\rm d}x= 1$. The first step gives rise to an entropy term that scales, for large N, as $O(N)$ (see for instance [30]). But since the energy $\mathcal{E}[\rho_w(x)]$ in equation (55) scales as N³, we can neglect the entropy term at leading order for large N. This gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z_N(w)\approx \int \mathcal{D}\rho_w ~{\rm exp}\left (- \beta \mathcal{E}[\rho_w(x)]\right)\,\delta \left(\int_{-\infty}^w {\rm d}x~\rho_w(x)-1 \right) \; \label{partition-2} \nonumber \end{align} \tag{ 56 }$$

where $\mathcal{D}\rho_w$ denotes the measure of a functional integral over all possible densities satisfying the normalisation constraint $\int_{-\infty}^w {\rm d}x~\rho_w(x) = 1$. To proceed further we replace the delta function by its integral representation and get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle Z_N(w)=N^3 \int \frac{{\rm d} \mu}{2 \pi {\rm i}} \int \mathcal{D}\rho_w ~{\rm exp}\left (- N^3 \,\mathcal{S}[\rho_w(x)] \right), ~{\rm with}\label{partition-3} \nonumber \end{align} \tag{ 57 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \mathcal{S}[\rho_w(x)]=& \left[\frac{1}{2} \int_{-\infty}^w {\rm d}x~x^2\rho_w(x) - \al \int_{-\infty}^w {\rm d}x\int_{-\infty}^w {\rm d}y~\rho_w(x)\rho_w(y)~|x-y| \right. \nonumber \\ &\left. +~\mu \left(\int_{-\infty}^w {\rm d}x~\rho_w(x)-1 \right) \right]. \label{action} \nonumber \end{align} \tag{ 58 }$$

The integral in (57) can be performed, for large N, by a saddle point approximation that gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Z_saddle} Z_N(w) \approx {\rm exp}\left (- N^3 S[\rho^*_w(x)] \right), \nonumber \end{align} \tag{ 59 }$$

where $\rho^*_w(x)$ is the saddle point density that minimises the action $S[\rho_w(x)]$ in (58). The equation for $\rho^*_w(x)$ is obtained from $\newcommand{\f}{f} \left(\frac{\delta S[\rho_w]}{\delta \rho_w(x)}\right)_{\rho_w=\rho^*_w}=0$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \frac{1}{2}x^2 -2\alpha \int_{-\infty}^w {\rm d}y~\rho^*_w(y)~|x-y| +\mu = 0 . \label{sadl-eq-1} \nonumber \end{align} \tag{ 60 }$$

This equation holds for $-\infty \leqslant x \leqslant w$. But, the support of $\rho_w^*(x)$ can not extend over the full range $(-\infty, w]$ due to the following reasons: when $x \to -\infty$, the first term in (60) grows as x², while the second term grows as $|x|$—hence they can not compensate each other. However the equation (60) holds true. The only possible way this can happen is that $\rho^*_w(x)$ has a finite support, say over $[-B, w]$ where B can be determined from the normalisation constraint

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{-B}^{w} \rho^*_w(x) {\rm d}x =1. \label{norml-1} \nonumber \end{align} \tag{ 61 }$$

Outside this region $\rho_w^*(x)$ is zero. For $x \in [-B, w]$, differentiating twice the saddle point equation (60) and using the identity $\newcommand{\f}{f} \frac{{\rm d}^2}{{\rm d}x^2}|x-y|=2\delta(x-y)$, it is easy to show that $\newcommand{\al}{\alpha} \rho^*_w(x) = 1/(4 \alpha)$. If $\newcommand{\al}{\alpha} w > 2 \al$, the saddle point density is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{rho_star_eq} \rho_w^*(x) = \frac{1}{4 \al} \; \quad {\rm for} \; \quad -2\al \leqslant x \leqslant 2 \al , \quad \quad w > 2\al. \nonumber \end{align} \tag{ 62 }$$

Thus for $\newcommand{\al}{\alpha} w > 2\al$ the charge density does not change from its flat equilibrium density—this is because the charges do not feel the presence of the wall. However, when $\newcommand{\al}{\alpha} w<2\al$, the wall tries to push the charges to the left of $\newcommand{\al}{\alpha} 2\al$ (see the left panel of figure 3). We have seen from above that the bulk density does not change from its equilibrium value $\newcommand{\al}{\alpha} \rho^*_w (x) = 1/(4\al)$ to the left of the wall at w. Normalisation to unity of the charge density then implies that the extra charge that the wall displaces must be accumulated at the wall, since the bulk is not affected. This leads, for $\newcommand{\al}{\alpha} w<2\al$, to a new saddle point density of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \rho^*_w(x)= \frac{1}{4\al}+ C~\delta(x-w),~{\rm for}~-B\leqslant x < w, \label{rho-guess} \nonumber \end{align} \tag{ 63 }$$

where C represents the density of the charges displaced and absorbed at the wall. We have three unknowns: $B, ~C$ and μ which are to be determined now. From the normalisation condition $\int_{-B}^w \rho^*_w(x)\, {\rm d}x = 1$ we get the relation between the two parameters B and C via

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \frac{(w+B)}{4\al} + C=1 . \label{eq-3} \nonumber \end{align} \tag{ 64 }$$

Zoom In Zoom Out Reset image size

We need two more equations. For that, we substitute the saddle point density $\rho^*_w(x)$ from (63) in (60) to get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \frac{1}{2}x^2 -\frac{1}{2} \int_{-B}^w {\rm d}y~|x-y|-2 \al~C~ |x-w| +\mu = 0 . \label{eq-17-supp} \nonumber \end{align} \tag{ 65 }$$

Now performing the integral over y explicitly, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \left(2\al~C + \frac{w-B}{2}\right)x ~+~\left(\mu - 2 \al~C~w-\frac{B^2+w^2}{4}\right)=0, ~{\rm for}~-B\leqslant x <w . \label{eq-18-supp} \nonumber \end{align} \tag{ 66 }$$

Since (66) is valid for arbitrary $x \in [-B, w]$, the coefficients of different powers of x are individually zero. As a result we get two additional equations

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle 2\al\,C + \frac{w-B}{2}=0, \label{eq-1} \nonumber \end{align} \tag{ 67 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \mu = 2 \al~C~w+\frac{B^2+w^2}{4}. \label{eq-2} \nonumber \end{align} \tag{ 68 }$$

We therefore have three equations (64), (67) and (68) for three unknowns $\mu, B$ and C, solving which we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle B=2\al, \label{B} \nonumber \end{align} \tag{ 69 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle C=\frac{1}{2} - \frac{w}{4\al}, \label{A} \nonumber \end{align} \tag{ 70 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \mu = \al^2 + \al w - \frac{w^2}{4} . \label{mu} \nonumber \end{align} \tag{ 71 }$$

Since $B, ~w$ and α are all positive by definition, then normalisation condition in (64) implies $C \leqslant 1$. This indicates that the above analysis is valid only for $\newcommand{\al}{\alpha} w > -2\alpha$. When $\newcommand{\al}{\alpha} w \to -2\al$, $C \to 1$: this means that all the charges are absorbed at the wall and there is no bulk charge left. Thus for $\newcommand{\al}{\alpha} w < -2\al$, we have effectively a single charge located at w subjected to a harmonic potential. Therefore, for the saddle point density $\rho_w^*(x)$ we have the following expressions, valid for all w

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{rho_summary} \rho^*_w(x) = \left\{\begin{array}{@{}lclc@{}} \frac{1}{4\al}, &-2\al \leqslant x \leqslant 2\al , & {\rm for} & w > 2\al \nonumber \\ \frac{1}{4\al} + \left(\frac{1}{2} - \frac{w}{4\al} \right)~\delta(x-w), & -2\al\leqslant x \leqslant w , & {\rm for} & -2\al \leqslant w \leqslant 2\al \nonumber \\ \delta(x-w) & & {\rm for} & w < -2\al . \end{array}\right. \nonumber \end{align} \tag{ 72 }$$

This saddle point density $\rho^*_w(x)$ has a nice interpretation. When the wall position $\newcommand{\al}{\alpha} w > 2\al$, it is given by the unperturbed density given in the first line of equation (72)—the charges do not feel the presence of the wall. When the wall position $\newcommand{\al}{\alpha} -2\al<w<2\al$, the wall displaces the charges over the region $\newcommand{\al}{\alpha} [w, 2 \al]$ and absorbs them on the wall as shown by the delta function term in the second line of equation (72) (see also the left panel of figure 3). Finally, when $\newcommand{\al}{\alpha} w<-2\al$, all the bulk charges are absorbed on the wall and the density is a simple delta function, given by the third line of equation (72).

Our next task is to insert this saddle point density in the action $\mathcal{S}$ in (58) and get the partition function Z_N(w) in (59) to leading order. Let us first consider $\newcommand{\al}{\alpha} w > 2 \al$. In this case $\newcommand{\al}{\alpha} \rho_w^*(x) = 1/(4 \al)$ for $\newcommand{\al}{\alpha} x \in [-2\al, +2\al]$. Substituting this density in (58) we get the saddle point action

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle S[\rho_w^*(x)] = - \frac{2}{3}\alpha^2, \quad \quad {\rm for} \quad w > 2\al . \label{S_largew} \nonumber \end{align} \tag{ 73 }$$

Therefore from (59) the partition function Z_N(w) for large N and for $\newcommand{\al}{\alpha} w > 2\al$ behaves as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle Z_N(w) \approx \ee^{\frac{2}{3}\alpha^2\,N^3} , \quad \quad {\rm for} \quad w > 2\al . \label{Z_largew} \nonumber \end{align} \tag{ 74 }$$

In particular, taking $w \to \infty$ limit, we obtain the denominator in (34) as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{denominator} Z_N(\infty) \approx \ee^{\frac{2}{3}\alpha^2\,N^3} . \nonumber \end{align} \tag{ 75 }$$

Hence, finally, for $\newcommand{\al}{\alpha} w > 2\al$, to leading order for large N, we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{Q_largew} Q(w,N) = \frac{Z_N(w)}{Z_N(\infty)} \approx 1 , \quad \quad {\rm for} \quad w > 2\al . \nonumber \end{align} \tag{ 76 }$$

To calculate the corrections to this leading order result, we need to consider the right large deviations function, that will be computed in the next section. Let us now consider the region where $\newcommand{\al}{\alpha} -2\al \leqslant w \leqslant 2\al$. Substituting the saddle point density $\rho_w^*(x)$ from the second line of (72) in (58) we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{S_inter} S[\rho_w^*(x)] = - \frac{8\al^3+12\al^2w-6\al w^2+w^3}{24 \al} , \quad \quad {\rm for} \quad -2\al \leqslant w \leqslant 2\al . \nonumber \end{align} \tag{ 77 }$$

Substituting this result in (59) and using the expression for the denominator in (75) we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{Q_mediumw} Q(w,N) = \frac{Z_N(w)}{Z_N(\infty)} \approx \ee^{- N^3\Phi_-(w)} , \nonumber \end{align} \tag{ 78 }$$

where the large deviation function $\Phi_-(w)$ actually has a very simple expression

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{Q_mediumw2} \Phi_-(w) = \frac{(2\al - w)^3}{24\,\al},~{\rm for} ~ -2\al \leqslant w \leqslant 2\al . \nonumber \end{align} \tag{ 79 }$$

Finally, we consider the region where $\newcommand{\al}{\alpha} w\leqslant-2\al$. In this case, substituting the saddle point density $\rho_w^*(x)$ from the third line of (72) in (58) we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{S_small} S[\rho_w^*(x)] = \frac{w^2}{2}, \quad \quad {\rm for} \quad w < -2\al . \nonumber \end{align} \tag{ 80 }$$

Substituting this result in (59) and using the expression for the denominator in (75) we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{Q_smallw} Q(w,N) &= \frac{Z_N(w)}{Z_N(\infty)} \approx \ee^{-N^3\Phi_-(w)} , \quad {\rm for}~\beta_0 \sim N \nonumber \\ \quad {\rm where,} & \Phi_-(w) = \frac{w^2}{2} + \frac{2}{3}\alpha^2,~{\rm for}~w\leqslant -2\al. \nonumber \end{align} \tag{ 81 }$$

In summary

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{sm-left_rf} \Phi_{-}(w) = \left\{\begin{array}{@{}cc@{}} \frac{(2\al-w)^3}{24\al},& -2\al\leqslant w \leqslant 2\al \nonumber \\ \frac{w^2}{2} + \frac{2}{3} \al^2 , & w \leqslant -2\al . \end{array}\right. \nonumber \end{align} \tag{ 82 }$$

In figure 4, we show a plot of $\Phi_-(w)$ as a function of w.

Zoom In Zoom Out Reset image size

Third order phase transition at $\newcommand{\al}{\alpha} w = 2\al$: the cumulative distribution $Q(w, N)$ in (34) is the ratio of two partition functions. Hence $-\ln Q(w, N) = -\ln Z_N(w)+\ln Z_N(\infty)$ can be interpreted as a free energy difference. Indeed, from (76), we see that to leading order for large N, $-\ln Q(w, N) \approx 0$ for $\newcommand{\al}{\alpha} w>2\al$. In contrast, for $\newcommand{\al}{\alpha} w<2\al$, using (78) and (81), we see that $-\ln Q(w, N) \approx N^3 \Phi_-(w)$ where $\Phi_-(w)$ is given in (82). Hence, we get (see figure 4)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{log_Qn} - \lim_{N \to \infty} \frac{\ln Q(w,N)}{N^3} = \left\{\begin{array}{@{}lcc@{}} 0 , & {\rm for}\; & w > 2\al \nonumber \\ \Phi_-(w) , & {\rm for}\;& w < 2\al . \end{array}\right. \nonumber \end{align} \tag{ 83 }$$

Thus $\Phi_-(w)$ is just the free energy cost in pushing the wall w to the left of the right edge $\newcommand{\al}{\alpha} 2 \al$ (see figure 3). From the expression of $\Phi_-(w)$ in the first line of (82), it follows that $\Phi_-(w)$ vanishes as the third power $\newcommand{\al}{\alpha} \Phi_-(w) \propto (2\al-w){}^3$ as $\newcommand{\al}{\alpha} w \to 2\al$ from the left. Thus the third derivative of the free energy is discontinuous at the critical point $\newcommand{\al}{\alpha} 2 \al$, making this a third order phase transition. Indeed, the pressure on the wall $P = -N^3 \Phi'_-(w)$ (derivative of the free energy with respect to the wall position) is zero for $\newcommand{\al}{\alpha} w > 2\al$ (the charges do not touch the wall) and is non zero for $\newcommand{\al}{\alpha} -2\al<w<2\al$. The mechanism of this third order transition is thus similar to the log-gas case [33]. However, in contrast to the log gas case, there is an additional third-order phase transition in the jellium model when $\newcommand{\al}{\alpha} w \to -2\al$ (see equation (82)). Indeed the third derivative of $\Phi_-(w)$ in equation (82) is also discontinuous at $\newcommand{\al}{\alpha} w = - 2 \al$. This transition is not of the 'pushed–pulled' type like the one at $\newcommand{\al}{\alpha} w =2\al$, but rather a condensation-type transition as all charges accumulate at the wall for $\newcommand{\al}{\alpha} w \leqslant -2\al$.

Interestingly, a similar third-order phase transition between the pushed and the pulled phase was recently found [40] by analysing large deviation functions associated with the position of the farthest charge in a d-dimensional jellium model. The limiting distribution of the position of the farthest charge is known in d  =  1 (and was computed by Baxter, see equations (32) and (47)) and in d  =  2 where the distribution, properly centred and scaled, approaches a Gumbel distribution [55]. However, for d  >  3, no explicit result is known for this limiting distribution. In d  =  1 this corresponds to the distribution of the maximum of |x_i|'s of the charges, as discussed above (see equations (32) and (48)). It was further shown that this observable exhibits a similar third order phase-transition even for short-range interactions, like the Yukawa potential [41], in $\newcommand{\g}{g} d \geqslant 1$.

4.2.2. Right large deviation.

We now focus on the distribution $Q(w, N)$, for large N, in the region $\newcommand{\al}{\alpha} 0<w-2\al \sim O(1)$, that characterises the large fluctuations of order $O(1)$ to the right of the mean. From the analysis performed in the previous section, we have seen that in this regime, to leading order for large N, $Q(w, N) \approx 1$ (see (76)). To compute the sub leading corrections to this leading order term 1, it is convenient to consider instead the PDF of $x_{\max}$, given by the derivative of (34)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{exact_PDF} P(w,N) =&~ \partial_w Q(w,N) = \frac{N}{Z_N(\infty)} \ee^{-\frac{N^2}{2}w^2} \int_{-\infty}^w {\rm d}x_1 \cdots \int_{-\infty}^w {\rm d}x_{N-1}\; \nonumber \\ &\times \exp \left(2\alpha\,N\, \sum_{j=1}^{N-1} |w-x_j| + \al\, N\,\sum_{1\leqslant i\neq j\leqslant N-1}|x_i-x_j| - \frac{N^2}{2}\sum_{i=1}^{N-1}x_i^2 \right) \nonumber \end{align} \tag{ 84 }$$

where we have simply separated out the x_N  =  w from the rest in equations (35) and (13). This can be re-written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{PDF_largeN_1} P(w,N) = \frac{N\,Z_{N-1}(\infty)}{Z_N(\infty)}\, \ee^{-\frac{N^2}{2}w^2} \Big \langle \ee^{2\al \,N\sum\nolimits_{j=1}^{N-1}(w-x_j)}\Big\rangle_{N-1} , \nonumber \end{align} \tag{ 85 }$$

where $\langle \cdots \rangle_{N-1}$ denotes the average over the Boltzmann distribution of N  −  1 charges. We can then analyse this average for large N, for $\newcommand{\al}{\alpha} w > 2 \al$, following [31] for the log-gas in the corresponding right large deviation regime. To evaluate this average, we note that essentially one single charge out of N is detached at $\newcommand{\al}{\alpha} w > 2\al$, while the rest of N  −  1 charges should be in their equilibrium flat configuration, i.e. with a density $\newcommand{\al}{\alpha} \rho^*_w(x) = 1/(4 \al)$ for $\newcommand{\al}{\alpha} x \in [-2\al, 2\al]$ (see the right panel of figure 3). Furthermore, for large N, to leading order, we can (i) approximate the average of the exponential in equation (85) by the exponential of the average and (ii) use that, to leading order for large N, $\newcommand{\f}{f} \frac{N\, Z_{N-1}(\infty)}{Z_N(\infty)} \sim {\rm e}^{-C_0\, N^2}$ for some constant C₀ (independent of w) to write

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle \label{PDF_largeN_2} P(w,N) \approx \ee^{-\frac{N^2}{2}w^2 + 2 \al \, N \langle \sum\nolimits_{j=1}^{N-1} (w-x_j)\rangle-N^2\,C_0} \approx \ee^{-N^2\left(\frac{w^2}{2}-2\al \int_{-2\al}^{2\al} (w-x)\rho^*_w(x) {\rm d}x - C_0\right)} . \nonumber \end{align} \tag{ 86 }$$

Using $\newcommand{\al}{\alpha} \rho^*_w(x) = 1/(4 \al)$ for $\newcommand{\al}{\alpha} x \in [-2\al, 2\al]$ and performing the integral in (86), we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ee}{{\rm e}} \displaystyle \label{sm_Epulled} P(w,N) \approx \ee^{-\Delta E_{\rm pulled}} \approx \ee^{-N^2\, \Phi_+(w)} , \nonumber \end{align} \tag{ 87 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{sm-phi+} \Phi_+(w) = \frac{(w-2\al)^2}{2} , \quad w > 2\al . \nonumber \end{align} \tag{ 88 }$$

Thus $\Delta E_{\rm pulled}$ in (87) corresponds to the energy in pulling out a single charge from the equilibrium configuration of charges with a flat density.

In this section we study the PDF P_G(g, N) of the gap $g=x_N-x_{N-1}$ between the positions of the two rightmost charges with ordered positions x_N and $x_{N-1}<x_N$. We show that the typical fluctuations of g, of order $O(1/N)$ for large N, are described by the scaling form,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle P_G(g,N) &\approx N~h_\alpha(gN),~{\rm with}, \nonumber \\ h_\alpha(z) &= \Theta(z)\,A(\alpha)\,\int_{-\infty}^\infty {\rm d}u~{\rm e}^{-u^2/2}~\frac{{\rm d}}{{\rm d}u}F_\alpha(u-z+4\alpha) \nonumber \end{align} \tag{ 89 }$$

where $\newcommand{\al}{\alpha} F_\alpha(y)$ is given by the solution of (47) with its associated eigenvalue $\newcommand{\al}{\alpha} A(\alpha)$. For much larger values of g, i.e. $g = O(1)$, P_G(g, N) is described by the large deviation form given in equation (103).

The distribution of the gap $g=x_N-x_{N-1}$ can be formally written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle P_G(g,N)= \frac{N!}{Z_N} \int {\rm d}x_1\ldots {\rm d}x_N {\rm e}^{-\beta E[\{x_i\}]} ~\delta(x_N-x_{N-1}-g)~\prod_{j=2}^N \Theta(x_j - x_{j-1}) , \label{gap_dist-1} \nonumber \end{align} \tag{ 90 }$$

where $\beta E[\{x_i\}]$ is given in equation (27). In the ordered sector, rewriting the energy $E[\{x_i\}]$ as in (29) and performing the change of variables $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} x_i \to \ep_i$ as in (37), one can write

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle P_G(g,N)=N^{1-N}\frac{N!}{Z_N} \int_{-\infty}^{\infty} \prod_{k=1}^N {\rm d}\ep_k\, {\rm e}^{-\frac{1}{2}\sum\nolimits_{k=1}^N \ep_k^2}\, \left[\,\prod_{k=2}^N \Theta(\ep_k-\ep_{k-1}+4\alpha)\, \delta(\ep_N-\ep_{N-1}-gN+4\alpha) \right]. \nonumber \end{align} \tag{ 91 }$$

It is useful to regroup the pair of variables $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \ep_N$ and $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \ep_{N-1}$, keeping the rest N  −  2 variables together and rewrite the integral as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle &P_G(g,N)= N^{1-N}\frac{N!}{Z_N}\int_{-\infty}^{\infty} {\rm d}\ep_N \int_{-\infty}^{\infty} {\rm d}\ep_{N-1}{\rm e}^{-\frac{1}{2} (\ep_N^2+\ep_{N-1}^2)} \delta(\ep_N-\ep_{N-1}-gN+4\alpha)\, \nonumber \\ &~~~~~~~~~~~~~~\times~\Theta(\ep_N-\ep_{N-1}+4\alpha)~\left(\int_{-\infty}^{\infty} \prod_{k=1}^{N-2} {\rm d}\ep_k\, {\rm e}^{-\frac{1}{2}\sum\nolimits_{k=1}^{N-2} \ep_k^2}\, \prod_{k=2}^{N-1} \Theta(\ep_k-\ep_{k-1}+4\alpha)\right)\, \nonumber \\ &~= \frac{N!N^{1-N}}{Z_N}\int_{-\infty}^{\infty} {\rm d}\ep_N \int_{-\infty}^{\infty} {\rm d}\ep_{N-1}{\rm e}^{-\frac{1}{2} (\ep_N^2+\ep_{N-1}^2)} \delta(\ep_N-\ep_{N-1}-gN+4\alpha)\, \nonumber \\ &~~~~~~~~~~~~~\times~\Theta(\ep_N-\ep_{N-1}+4\alpha)~D_\alpha(\ep_{N-1}+4\alpha,N-2), \label{eq:PNg} \nonumber \end{align} \tag{ 92 }$$

where in the last step we have used the definition of $\newcommand{\al}{\alpha} D_\alpha(x, N)$ in (41). This formula for the PDF of the gap P_G(g,N) in equation (92) is exact for any N. We now analyze it in the large N limit. In this limit, it turns out that the typical fluctuations of the gap are of order $O(1/N)$, as suggested by the appearance of the scaling variable $g\, N$ in equation (92), while, as for $x_{\max}$, the atypical fluctuations are of order $O(1)$. We now analyse separately these two regimes of typical and atypical fluctuations of the first gap.

5.1. Typical fluctuations of the gap g

To analyse the typical fluctuations of the gap g, we consider the limit $N \to \infty$, $g \to 0$, keeping the scaling variable $z = N\, g$ fixed. For large N we use equation (42) to write $\newcommand{\al}{\alpha} D_\alpha(x, N-2)= D_\alpha(\infty, N-2) F_\alpha(x, N-2)$. Using further $\newcommand{\al}{\alpha} Z_N \approx N!D_\alpha(\infty, N)/N^N$, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle P_G(g,N) \approx &~ N\frac{D_\alpha(\infty,N-2)}{D_\alpha(\infty,N)}\int_{-\infty}^{\infty} {\rm d}\ep_N \int_{-\infty}^{\infty} {\rm d}\ep_{N-1}{\rm e}^{-\frac{1}{2} (\ep_N^2+\ep_{N-1}^2)} \delta(\ep_N-\ep_{N-1}-gN+4\alpha) \nonumber \\ &~~\times~\Theta(\ep_N-\ep_{N-1}+4\alpha)\,~F_\alpha(\ep_{N-1}+4\alpha,N-2). \nonumber \end{align} \tag{ 93 }$$

In the large N limit, we then use the fact that $\newcommand{\al}{\alpha} D(\alpha, N)\sim [A(\alpha)]^{-N}$ and that $\newcommand{\al}{\alpha} F_\alpha(x, N \to \infty) = F_\alpha(x)$. Further, keeping $N\, g =z$ fixed in the scaling limit, we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle P_G(g,N) \approx NA(\alpha)^2\int_{-\infty}^{\infty} {\rm d}x \int_{-\infty}^{\infty} {\rm d}y~{\rm e}^{-\frac{1}{2} (x^2+y^2)} \delta(x-y-gN+4\alpha) ~\Theta(x-y+4\alpha)\,~F_\alpha(y+4\alpha), \label{PN_g} \nonumber \end{align} \tag{ 94 }$$

where $\newcommand{\al}{\alpha} F_\alpha(x)$ satisfies the differential equation (47). The expression for P_G(g,N) in equation (94) can be written in the scaling form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle \label{scaling_PNg1} P_G(g,N) \approx N\,h_\alpha(g\,N) , \nonumber \end{align} \tag{ 95 }$$

where the scaling function $\newcommand{\al}{\alpha} h_\alpha(z)$ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{h_a_1} h_\alpha(z) = A(\alpha)^2\int_{-\infty}^{\infty} {\rm d}x \int_{-\infty}^{\infty} {\rm d}y~{\rm e}^{-\frac{1}{2} (x^2+y^2)} \delta(x-y-z+4\alpha) ~\Theta(x-y+4\alpha)\,~F_\alpha(y+4\alpha). \nonumber \end{align} \tag{ 96 }$$

This scaling function $\newcommand{\al}{\alpha} h_\alpha(z)$ in equation (96) can be further simplified as follows. The presence of the theta function in equation (96) indicates that this integral is non-zero only when $\newcommand{\al}{\alpha} x > y - 4 \alpha$. In contrast, the delta function indicates that this is non-zero only when $\newcommand{\al}{\alpha} x = y+z-4\al$. Hence, for z  <  0, these two conditions can not be satisfied simultaneously. This indicates that $\newcommand{\al}{\alpha} h_\alpha(z<0) = 0$. For z  >  0, once the delta function constraint is satisfied, then the theta function constraint is automatically satisfied. Hence, for z  >  0, we can write

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{h_a_2} h_\alpha(z) = \Theta(z)\,A^2(\alpha)\, \int_{-\infty}^{\infty} {\rm d}x \int_{-\infty}^{\infty} {\rm d}y~{\rm e}^{-\frac{1}{2} (x^2+y^2)} \delta(x-y-z+4\alpha)~F_\alpha(y+4\alpha) . \nonumber \end{align} \tag{ 97 }$$

Performing the integral over x, we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{h_a_3} h_\alpha(z) = A^2(\alpha)\, \Theta(z) \int_{-\infty}^\infty {\rm d}y\, {\rm e}^{-\frac{y^2}{2}} {\rm e}^{-\frac{(y+z-4\al)^2}{2}}~F_\alpha(y+4\alpha) . \nonumber \end{align} \tag{ 98 }$$

Using the differential equation (47) one can simplify further

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{h_a_4} h_\alpha(z) = A(\alpha) \Theta(z) \int_{-\infty}^\infty {\rm d}y\, {\rm e}^{-\frac{1}{2}(y+z-4\al)^2} F_\alpha'(y) . \nonumber \end{align} \tag{ 99 }$$

One can also do an integration by parts to rewrite it as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{h_a_5} h_\alpha(z) = A(\alpha) \Theta(z) \int_{-\infty}^\infty {\rm d}y \, (y+z-4\al) {\rm e}^{-\frac{1}{2}(y+z-4\al)^2} F_\alpha(y) . \nonumber \end{align} \tag{ 100 }$$

In figure 5 we compare this theoretical result with numerical simulation and observe a very good agreement. One can also estimate the asymptotic tails of the scaling function $\newcommand{\al}{\alpha} h_\alpha(z)$. From equation (100), as $z \to 0$, the scaling function $\newcommand{\al}{\alpha} h_\alpha(z)$ approaches a constant given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{h0} h_\alpha(0) = A_\alpha \int_{-\infty}^\infty {\rm d}y \, (y-4\al) {\rm e}^{-\frac{1}{2}(y-4\al)^2} F_\alpha(y) , \nonumber \end{align} \tag{ 101 }$$

which can be evaluated numerically. For $z \to \infty$, one can show, by analysing the integral in equation (100) and using the tails of $\newcommand{\al}{\alpha} F_\alpha(y)$ from equation (52), that to leading order for large z,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle \label{h_largez} h_\alpha(z) \sim {\rm e}^{-z^2/2 + o(z^2)} . \nonumber \end{align} \tag{ 102 }$$

Zoom In Zoom Out Reset image size

As expected, this is similar to the right tail behaviour of the PDF of $x_N = x_{\max}$ (see equation (52)), since to create a large gap, we must have $\newcommand{\g}{g} x_N \gg x_{N-1}$. The fact that the large gap asymptotic behaviour coincides with the right tail of $x_{\max}$ also holds for the log-gas case [44].

5.2. Atypical large fluctuations of the gap g

To analyse the large fluctuations of the gap g of order $O(1)$, it is useful to remark that the configurations that contribute to the PDF $P(g, N)$ are the same that contribute to a large value of $x_{\max}$ to the right of its mean, as depicted in the right panel of figure 3. In such configurations, $x_N = x_{\max} = w$ while the second particle x_N−1 is located close to the right edge $\newcommand{\al}{\alpha} x_{N-1} \approx 2 \alpha$, leading to a gap $\newcommand{\al}{\alpha} g = w-2\alpha$. Therefore, for large N, one obtains that $\newcommand{\al}{\alpha} {\rm Prob.}(x_{N}-x_{N-1} = g) \approx {\rm Prob.}(x_{\max} = g + 2 \al)$. Therefore, from the right large deviation form of the PDF of $x_{\max}$ obtained in equations (87) and (88), one gets, the large deviation form of $P(g, N)$ to leading order for large N as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{large_dev_gap} P(g,N) \sim {\rm e}^{-N^2 \Psi_+(g)} , \; \Psi_+(g) = \Phi_+(g + 2\al) = \frac{g^2}{2} , \nonumber \end{align} \tag{ 103 }$$

where we have used the explicit expression of $\Phi_+(w)$ given in equation (88). Thus, this large deviation regime matches exactly with the right tail of the regime of typical fluctuations (see equations (95) and (102)).

In this section, we study the statistics of the index $N_+ = \sum\nolimits_{i=1}^N \Theta(x_i)$. This random variable lies within the range $0\leqslant N_+\leqslant N$ and we now compute its distribution $P_I(N_+, N)$ for large N. It is also clear that $\langle N_+\rangle=N/2$ and the distribution $P_I(N_+, N)$ is symmetric around this mean. Given the joint PDF $\mathcal{P}(x_1, x_2, \ldots, x_N)$ in (12) along with (14), $P_I(N_+, N)$ can be expressed as a multiple integral

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_I(N_+,N)= \int_{-\infty}^{\infty}{\rm d}x_1\dots \int_{-\infty}^{\infty}{\rm d}x_N\, \mathcal{P}(x_1,x_2,\dots,x_N)\, \delta\left[\sum_{i=1}^N \Theta(x_i)-N_+\right]. \label{index_dist.1} \nonumber \end{align} \tag{ 104 }$$

Since the integrand in (104) is symmetric under any permutation of the x_i's, we can order the x_i's, with $x_1<x_2< x_3\dots < x_N$ and rewrite it as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle P_I(N+,N)= \frac{N!}{Z_N} \int {\rm d}x_1\ldots {\rm d}x_N \ee^{-\beta\,E[\{x_i\}]} \delta\left[\sum_{i=1}^N \Theta(x_i)-N_+\right]\, \prod_{j=2}^N \Theta(x_j - x_{j-1}) \label{index_dist.2} , \nonumber \end{align} \tag{ 105 }$$

where $\beta\, E[\{x_i\}]$ is given in equation (14). In this ordered sector $(x_1<x_2<x_3\ldots<x_N)$, we use the same trick to eliminate the absolute values as done in (29). Hence, up to an overall normalisation constant, we can write

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle P_I(N_+,N) \propto & \int_{-\infty}^{\infty}{\rm d}x_1\dots \int_{-\infty}^{\infty} {\rm d}x_N \, {\rm e}^{-\frac{N^2}{2}\sum\nolimits_{i=1}^N [x_k- \frac{2\alpha}{N} (2k-N-1)]^2}\, \nonumber \\ &~~~~\times~\prod_{j=2}^N \Theta(x_j - x_{j-1}) \, \delta\left[\sum_{i=1}^N \Theta(x_i)-N_+\right]\, . \label{index_dist.3} \nonumber \end{align} \tag{ 106 }$$

Next we perform the change of variables given in (37) and rewrite the product of theta functions in equation (106) as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \displaystyle \prod_{k=2}^N \Theta(x_k-x_{k-1})= \prod_{k=2}^N \Theta(\ep_k-\ep_{k-1}+4\alpha)\, . \label{constraint.1} \nonumber \end{align} \tag{ 107 }$$

We want to compute the probability that N₊ charges are on the positive side, or equivalently that $N_{-}=N-N_+$ charges are on the negative side. Consider first the N₋ charges on the negative side. In the ordered sector, we have to ensure that the position of the N₋'th charge is negative. This automatically ensures (since we are in the ordered sector) that all the N₋ charges with positions x₁, x₂, ..., $x_{N_{-}}$ are negative. Thus, using the variables $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \epsilon_i$'s in equation (37), this condition translates to (see figure 6)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle x_{N_{-}}<0 \quad \quad {\rm or~~equivalently}\quad \quad \ep_{N_{-}}< 2\alpha(N+1-2N_{-})= 4\alpha \left(N_+-\frac{N}{2}\right)+ 2\alpha \label{cond.1} \nonumber \end{align} \tag{ 108 }$$

where we have used (37) with i  =  N₋ and $N_{-}=N-N_+$. Similarly, the condition $x_{N_{-}+1}>0$ automatically ensures (in the ordered sector) that the position of all N₊ charges on the right are positive, i.e. $x_{N_{-}+1}>0$, $x_{N_{-}+2}>0$, ..., x_N  >  0. Thus this condition translates to (see figure 6)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle x_{N_{-}+1}>0 \quad \quad {\rm or~~equivalently}\quad \quad \ep_{N_{-}+1}> 2\alpha(N-1-2N_{-}) = 4\alpha \left(N_+-\frac{N}{2}\right)- 2\alpha\, . \label{cond.2} \nonumber \end{align} \tag{ 109 }$$

Zoom In Zoom Out Reset image size

For later convenience, let us define

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle z= 4\alpha \left(N_+- \frac{N}{2}\right)\, . \label{def_z} \nonumber \end{align} \tag{ 110 }$$

Then, in terms of the z variable, the two conditions in (108) and (109) are expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \displaystyle \ep_{N_{-}}< z+ 2\alpha \quad \quad {\rm and}\quad \quad \ep_{N_{-}+1}> z- 2\alpha \, . \label{cond.3} \nonumber \end{align} \tag{ 111 }$$

Thus finally, using these new variables, $P_I(N_+, N)$ in (106) simplifies to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle &P_I(N_+,N) \propto T(N_+,N),~~{\rm where},\nonumber \\ &T(N_+,N)= \int_{-\infty}^{\infty} \prod_{k=1}^N {\rm d}\ep_k\, \ee^{-\frac{1}{2}\sum\nolimits_{k=1}^N \ep_k^2}\, \prod_{k=2}^N \Theta(\ep_k-\ep_{k-1}+4\alpha)\, \nonumber \\ &~~~~~~~~~~~~~~~~~~\times~\Theta(z+2\alpha- \ep_{N_{-}})\, \Theta(\ep_{N_{-}+1}- z+2\alpha)\nonumber , \label{index_dist.4} \nonumber \end{align} \tag{ 112 }$$

where we have used the condition in (107), as well as the two conditions in (111). Note that these equations are strictly valid for 0  <  N₊  <  N. For N₊  =  N and N₊  =  0, one has to consider a slightly different integral, but this does not make any difference in the scaling limit. Thus, once again we have reduced our original problem of a long-ranged Coulomb gas to a problem of a short-ranged gas where there are only nearest neighbour interactions. Additionally, now there is a 'defect' on the bond connecting $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \ep_{N_{-}}$ and its right neighbour $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \ep_{N_{-}+1}$ that makes this short-ranged gas inhomogeneous (see figure 6).

The integral in (112) can be further simplified into two blocks as follows. For simplicity, let us denote the value of $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \ep$'s across the 'defect' as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \displaystyle \ep_{N_{-}}=u \quad \quad {\rm and} \quad \quad \ep_{N_{-}+1}=v\, . \label{def_uv} \nonumber \end{align} \tag{ 113 }$$

Then the integral in (112) can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle T(N_+,N) = \int_{-\infty}^{\infty} {\rm d}u \int_{-\infty}^{\infty} {\rm d}v \, T_1(u)\, T_2(v)\, \ee^{-\frac{1}{2} (u^2+v^2)}\, \Theta(v-u+4\alpha)\, \Theta(z+2\alpha- u)\, \Theta(v-z+2\alpha)\, , \label{index_dist.5} \nonumber \end{align} \tag{ 114 }$$

where T₁(u) is the integral over the left M  =  N₋  −  1 variables (given $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \ep_{N_{-}}=u$) and $T_2(v)$ is the integral over the right N₊  −  1  =  N  −  M  −  2 variables (given $\newcommand{\ep}{\epsilon} \newcommand{\e}{{\rm e}} \ep_{N_{-}+1}=v$). They are given explicitly by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle T_1(u)= \int_{-\infty}^{\infty}{\rm d}\ep_1\int_{-\infty}^{\infty} {\rm d}\ep_2\dots \int_{-\infty}^{\infty} {\rm d}\ep_{M}\, \ee^{-\frac{1}{2} \sum\nolimits_{k=1}^{M} \ep_k^2}\, \left[\prod_{k=2}^{M}\Theta(\ep_k-\ep_{k-1}+4\alpha)\right]\, \Theta(u-\ep_M+4\alpha)\, \label{T1u} \nonumber \end{align} \tag{ 115 }$$

where we recall that M  =  N₋  −  1. Similarly, for the right block, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle T_2(v)= & \int_{-\infty}^{\infty}{\rm d}\ep_{M+3}\int_{-\infty}^{\infty} {\rm d}\ep_{M+4}\dots \int_{-\infty}^{\infty} {\rm d}\ep_{N}\,~\ee^{-\frac{1}{2} \sum\nolimits_{k=M+3}^{N} \ep_k^2}\, \nonumber \\ & \times ~\left[\prod_{k=M+4}^{N}\Theta(\ep_k-\ep_{k-1}+4\alpha)\right]\, \Theta(\ep_{M+3}-v+4\alpha). \label{T2v} \nonumber \end{align} \tag{ 116 }$$

We can re-write T₁(u) in (115) by incorporating the constraints imposed by the theta functions directly in the limits of integration as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle T_1(u)= \int_{-\infty}^{u+4\alpha} {\rm d}\ep_M\, \ee^{-\frac{1}{2}\ep_M^2}\,\int_{-\infty}^{\ep_M+4\alpha} {\rm d}\ep_{M-1} \, \ee^{-\frac{1}{2}\ep_{M-1}^2}\ldots \int_{-\infty}^{\ep_2+4\alpha} {\rm d}\ep_1\, \ee^{-\frac{1}{2}\ep_{1}^2}. \label{T1u.1} \nonumber \end{align} \tag{ 117 }$$

Note that this is exactly, the function $\newcommand{\al}{\alpha} D_\alpha(x, N)$ defined in (41). Hence, T₁(u) in (117) reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle T_1(u)= D_\alpha(u+4\alpha, M)\, , \quad \quad {\rm where}\quad M= N_{-} - 1= N-N_+-1. \label{T1u.2} \nonumber \end{align} \tag{ 118 }$$

In a similar way, we can rewrite the integral $T_2(v)$ in (116) as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ep}{\epsilon} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle T_2(v)= \int_{v-4\alpha}^{\infty} {\rm d}\epsilon_{N_{-}+2}\, \ee^{-\frac{1}{2}\epsilon_{N_{-}+2}^2}\, \int_{\epsilon_{N_{-}+2}-4\alpha}^{\infty} {\rm d}\epsilon_{N_{-}+3}\, \ee^{-\frac{1}{2}\epsilon_{N_{-}+3}^2}\ldots \int_{\epsilon_{N-1}-4\alpha}^{\infty} {\rm d}\epsilon_N\, \ee^{-\frac{1}{2}\epsilon_N^2}\, . \label{T2v.1} \nonumber \end{align} \tag{ 119 }$$

As in the case of the left block, let us define a function similar to $\newcommand{\al}{\alpha} D_\alpha(x, M)$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle E_\alpha(x, M)= \int_{x}^{\infty} {\rm d}y_1\, \ee^{-\frac{1}{2} y_1^2}\, \int_{y_1-4\alpha}^{\infty} {\rm d}y_2\, {\rm e}^{-\frac{1}{2} y_2^2} \ldots \int_{y_{M-1}-4\alpha}^{\infty}{\rm d}y_M\, \ee^{-\frac{1}{2} y_M^2}\, . \label{Ex_def} \nonumber \end{align} \tag{ 120 }$$

Then, $T_2(v)$ in (119) can be written simply as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle T_2(v)= E_\alpha(v-4\alpha, N_+ -1)\, . \label{T2v.2} \nonumber \end{align} \tag{ 121 }$$

Finally, from the definitions of $\newcommand{\al}{\alpha} D_\alpha(x, M)$ (41) and $\newcommand{\al}{\alpha} E_\alpha(x, M)$ (120), it is easy to check, by performing the change of variables $y_k\to -y_k$, the following identity

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle D_\alpha(x,M)= E_\alpha(-x,M) \label{DE.0} \nonumber \end{align} \tag{ 122 }$$

valid for all $\newcommand{\g}{g} M\geqslant 0$. Plugging the results from equations (118) and (121) into (114) gives

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle T(N_+,N) =& \int_{-\infty}^{z+2\alpha} {\rm d}u\, \int_{z-2\alpha}^{\infty} {\rm d}v \,\, D_\alpha(u+4\alpha, N-N_+ -1)\,\, E_\alpha(v-4\alpha, N_+ -1)\,\, \nonumber \\ &\times \ee^{-\frac{1}{2} (u^2+v^2)}\, \,\Theta(v-u+4\alpha)\, \label{index_dist.6} \nonumber \end{align} \tag{ 123 }$$

where $\newcommand{\f}{f} \newcommand{\al}{\alpha} z= 4\alpha \left(N_+-\frac{N}{2}\right)\,$ from (110). The double integral in (123) can be further simplified by making the following observation. Let us look at the range of integration of u. The theta function $\newcommand{\al}{\alpha} \Theta(v-u+4\alpha)$ demands that $\newcommand{\al}{\alpha} u<v+4\alpha$. Hence we can eliminate the theta function and write it as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle T(N_+,N) = \int_{z-2\alpha}^{\infty} {\rm d}v \,\, \int_{-\infty}^{\min (v+4\alpha, z+2\alpha)} {\rm d}u\, D_\alpha(u+4\alpha, N-N_+ -1)\,\, E_\alpha(v-4\alpha, N_+ -1)\,\, \ee^{-\frac{1}{2} (u^2+v^2)}\, . \label{index_dist.7} \nonumber \end{align} \tag{ 124 }$$

However, the lower limit of the $v$ integration implies $\newcommand{\al}{\alpha} v>z-2\alpha$. This means $\newcommand{\al}{\alpha} v+4\alpha>z+2\alpha$. Hence, we necessarily have, $\newcommand{\al}{\alpha} {\min (v+4\alpha, z+2\alpha)}= z+2\alpha$. Thus we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle T(N_+,N) = \int_{z-2\alpha}^{\infty} {\rm d}v \,\, \int_{-\infty}^{z+2\alpha}{\rm d}u\, D_\alpha(u+4\alpha, N-N_+ -1)\,\, E_\alpha(v-4\alpha, N_+ -1)\,\, \ee^{-\frac{1}{2} (u^2+v^2)}\, . \label{index_dist.8} \nonumber \end{align} \tag{ 125 }$$

Making further the change of variable $v\to -v$, we can write it as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle T(N_+,N) = \int_{-\infty}^{-z+2\alpha} {\rm d}v \,\, \int_{-\infty}^{z+2\alpha} {\rm d}u\, D_\alpha(u+4\alpha, N-N_+ -1)\,\, E_\alpha(-v-4\alpha, N_+ -1)\,\, \ee^{-\frac{1}{2} (u^2+v^2)}\, . \label{index_dist.9} \nonumber \end{align} \tag{ 126 }$$

Finally, using $\newcommand{\al}{\alpha} E_\alpha(-x, M)= D_\alpha(x, M)$ from (122), we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle &P_I(N_+,N)\propto T(N_+,N),~~{\rm with}\nonumber \\ &T(N_+,N)=\left[\int_{-\infty}^{z+2\alpha} D_\alpha(u+4\alpha, N-N_+-1)\, \ee^{-\frac{1}{2} u^2}\, {\rm d}u\right] \nonumber \\ &\qquad\qquad\quad\,\,\times~\left[\int_{-\infty}^{-z+2\alpha} D_\alpha(v+4\alpha, N_+-1)\, \ee^{-\frac{1}{2} v^2}\, {\rm d}v\right]\label{index_dist.10} \nonumber \end{align} \tag{ 127 }$$

where we recall $\newcommand{\f}{f} \newcommand{\al}{\alpha} z= 4\alpha \left(N_+- \frac{N}{2}\right)$. Note that the rhs of (127) is manifestly a symmetric function of N₊ around N₊  =  N/2.

This result in (127) is exact for all 1  <  N₊  <  N. We now analyse it in the large N limit. For large N, it turns out that the typical fluctuations of N₊ around its mean N₊  =  N/2 are of order $O(1)$ (i.e. $z = O(1)$), while the atypical fluctuations are of order $O(N)$ (i.e. $z = O(N)$). Below, we analyse separately the probability distribution of typical and atypical fluctuations.

6.1. Typical fluctuations of N₊ 

It is convenient to define the fraction c  =  N₊ /N with $0\leqslant c\leqslant 1$. In the typical regime, where $\newcommand{\al}{\alpha} N_+ = N/2 + z/(4\al)$ with $z = O(1)$. This amounts to consider the scaling limit $c\to 1/2$, $N\to \infty$, while keeping the product $\newcommand{\al}{\alpha} z= 4\alpha N(c-1/2)$ finite, $O(1)$. To analyse (127) in this scaling limit ($N\to \infty$ with z fixed), we follow the method of section 3 and rewrite the function $\newcommand{\al}{\alpha} D_\alpha(x, M)$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle F_\alpha(x,M)= \frac{D_\alpha(x,M)}{D_\alpha(\infty,M)}. \label{Fx_def} \nonumber \end{align} \tag{ 128 }$$

In section 3, we have shown that in the large M limit, $\newcommand{\al}{\alpha} D_\alpha(\infty, M) \sim [A(\alpha)]^{-M}$ where $\newcommand{\al}{\alpha} A(\alpha)$ can be interpreted as the free energy associated to the short-ranged gas whose partition function is given in equation (41). Furthermore, for large M, the function $\newcommand{\al}{\alpha} F_\alpha(x, M)$ converges to a M independent limiting function $\newcommand{\al}{\alpha} F_\alpha(x)$ (as stated in (46)) where $\newcommand{\al}{\alpha} F_\alpha(x)$ satisfies the nonlocal eigenvalue equation (47). We then replace $\newcommand{\al}{\alpha} D_\alpha(x, M)= D_\alpha(\infty, M) F_\alpha(x, M)$ in (127), take the scaling limit using (44) and (46) and obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\ee}{{\rm e}} \newcommand{\f}{f} \displaystyle P_I(N_+,N) \propto \left[\int_{-\infty}^{z+2\alpha} A(\alpha)~F_\alpha(u+4\alpha)\, \ee^{-\frac{1}{2} u^2}\, {\rm d}u\right] \left[\int_{-\infty}^{-z+2\alpha} A(\alpha)~F_\alpha(v+4\alpha)\, \ee^{-\frac{1}{2} v^2}\, {\rm d}v\right], \label{index_dist.11} \nonumber \end{align} \tag{ 129 }$$

where we have absorbed the prefactor $\newcommand{\al}{\alpha} A(\al){}^{-N}$ in the proportionality constant. Furthermore, by using (47), the integrals over u and $v$ can be performed explicitly (using $\newcommand{\al}{\alpha} F_\alpha(x \to -\infty) = 0$, see equation (51)). This gives

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \displaystyle P_I(N_+,N) \propto F_\alpha(z+2\alpha)\, F_\alpha(-z+2\alpha)\, . \label{index_result.1} \nonumber \end{align} \tag{ 130 }$$

The proportionality constant can be fixed using the overall normalisation $\sum\nolimits_{N_+=0}^N P_I(N_+, N)=1$.

Summarising, the random variable N₊ typically fluctuates on a scale of $O(1)$ around its mean value N/2. We find that as $N\to \infty$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle N_+ \to \frac{N}{2} + \frac{1}{4\alpha}\, z \label{typical_scaling} \nonumber \end{align} \tag{ 131 }$$

where the random variable z has a limiting N-independent distribution $\newcommand{\al}{\alpha} f_\alpha(z)$. In other words, the distribution $P_I(N_+, N)$ converges to a limiting scaling form in the large N limit

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle P_I(N_+, N) \to 4\alpha\, f_\alpha\left(4 \alpha \left(N_+-\frac{N}{2}\right)\right)\, , \label{limiting.1} \nonumber \end{align} \tag{ 132 }$$

where the scaling function $\newcommand{\al}{\alpha} f_\alpha(z)$ is given from (130) as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle f_\alpha(z)= \frac{F_\alpha(z+2\alpha)\, F_\alpha(-z+2\alpha)}{\int_{-\infty}^{\infty} {\rm d}z\, F_\alpha(z+2\alpha)\, F_\alpha(-z+2\alpha)}. \label{limiting.2} \nonumber \end{align} \tag{ 133 }$$

The function $\newcommand{\al}{\alpha} f_\alpha(z)$ is manifestly symmetric around z  =  0. In figure 7 we compare this analytical result with numerical simulations and observe excellent agreement. The asymptotic behaviour of $\newcommand{\al}{\alpha} f_\alpha(z)$ for large z can be easily derived using the asymptotic decay of $\newcommand{\al}{\alpha} \newcommand{\e}{{\rm e}} F_\alpha(z\to -\infty) \sim \exp[- |z|^3/{24\alpha}]$ (see equation (52)) and the fact that $\newcommand{\al}{\alpha} F_\alpha(z\to \infty)=1$. Plugging these asymptotics in (133) gives

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle f_\alpha(z) \sim \exp[- \frac{1}{24\alpha}\,|z|^3]\quad {\rm as}\quad |z|\to \infty \, . \label{asymp.1} \nonumber \end{align} \tag{ 134 }$$

Zoom In Zoom Out Reset image size

Thus the limiting distribution $\newcommand{\al}{\alpha} f_\alpha(z)$ in equation (133) is obviously non-Gaussian. This is at variance with the log-gas where the typical fluctuations of the index are known to be Gaussian (see equation (10)). Furthermore, as we will see below, this tail behaviour from the central regime matches smoothly with the large deviation behavior of N₊.

6.2. Atypical large fluctuations of the index N₊ 

In this section we study large deviations regime of $P_I(N_+, N)$ in (104) where N₊  −  N/2  =  O(N) in the large N limit. Our starting point is the exact expression for $P_I(N_+, N)$ in (105) which we write as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle P_I(N_+,N)&=\frac{N!~I(N_+,N)}{Z_N} , \nonumber \\ I(N_+,N)&=\int {\rm d}x_1\ldots {\rm d}x_N {\rm e}^{-\beta\, E[\{x_i\}]} \delta\left[\sum_{i=1}^N \Theta(x_i)-N_+\right], \label{I_N+} \nonumber \end{align} \tag{ 135 }$$

with $\beta \, E[\{x_i\}]$ given in (14). Hence, I(N₊ , N) can be interpreted as the partition function of the 1d jelllium under the external constraint that there are exactly N₊ particles on the positive side. As the function I(N₊  =  cN, N) is symmetric around c  =  1/2, we assume, for convenience, that $0\leqslant c \leqslant 1/2$.

To proceed further, we follow the same Coulomb gas method as explained in section 4.2, to compute the partition function I(N₊  =  cN, N). We first replace the multiple integrals over x_i's in equation (135) by a functional integral over possible densities

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho_I(x,N)=N^{-1}\sum_i\delta(x-x_i) , \label{den-1} \nonumber \end{align} \tag{ 136 }$$

where the subscript I refers to 'index'. The density $\rho_I$ (i) is normalised and (ii) satisfies the constraint that N₊  =  cN charges are on the positive axis, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{constraint_rho_I} (i) \;\; \int_{-\infty}^\infty {\rm d}x \, \rho_I(x,N) = 1 , \;\;\; (ii) \int_{-\infty}^\infty {\rm d}x \, \Theta(x) \rho_I(x,N) = c . \nonumber \end{align} \tag{ 137 }$$

Therefore, the partition function $I(cN, N)$ reads, to leading order for large N

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I(cN,N) \propto \int \mathcal{D}[\rho_I]~{\rm e}^{-N^3\Sigma_c[\rho_I]} , \label{I_c} \nonumber \end{align} \tag{ 138 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \Sigma_c[\rho_I]=&~\frac{1}{2} \int_{-\infty}^{\infty} {\rm d}x~x^2\rho_I(x) - \alpha \int_{-\infty}^{\infty} {\rm d}x \int_{-\infty}^{\infty} {\rm d}y~\rho_I(x)\rho_I(y)~|x-y| \nonumber \\ &~~~+ A_1\left(\int_{-\infty}^{\infty}{\rm d}x~\rho_I(x) - 1\right) + A_2\left(\int_{-\infty}^{\infty}{\rm d}x~\Theta(x)\rho_I(x) - c\right), \label{Sigma-1} \nonumber \end{align} \tag{ 139 }$$

where A₁ and A₂ are Lagrange multipliers to enforce the constraints satisfied by $\rho_I(x)$ (137).

In the large N limit, the functional integral in equation (138) is dominated by the charge density $\rho_I^*$ that minimise $\Sigma_c[\rho_I]$. Numerical simulations indicate (see figure 8) that $\rho_I^*$ has two distinct supports and is of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho_I^*(x) = \left\{\begin{array}{@{}lcc@{}} \rho_1(x) + \lambda~\delta(x),&{\rm for}&-B\leqslant x \leqslant 0 \nonumber \\ \rho_2(x)&{\rm for}&a\leqslant x \leqslant b \end{array}\right., \label{den-anticipate} \nonumber \end{align} \tag{ 140 }$$

where the constants $\lambda >0, ~B>0, ~a>0$ and b  >  0 have to be determined. Inserting this form of density in the functional $\Sigma_c[\rho]$ in (139) we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{equ141} \Sigma_c[\rho^*_I]= &~\frac{1}{2} \int_{-B}^{0} {\rm d}x~x^2\rho_1(x)+\frac{1}{2} \int_{a}^{b} {\rm d}x~x^2\rho_2(x)\nonumber \\ &~~- \alpha \left(\int_{-B}^{0} {\rm d}x\int_{-B}^{0} {\rm d}y~\rho_1(x)\rho_1(y)|x-y| + 2 \lambda \int_{-B}^{0} {\rm d}x~|x|~\rho_1(x) \right. \nonumber \\ &~~ + 2\int_{-B}^{0} {\rm d}x\int_{a}^{b} {\rm d}y\rho_1(x)\rho_2(y)|x-y| +2 \lambda \int_{a}^{b} {\rm d}x~|x|~\rho_2(x) \nonumber \\ &~~\left.+ \int_{a}^{b} {\rm d}x\int_{a}^{b} {\rm d}y~\rho_2(x)\rho_2(y)|x-y| \right) \nonumber \\ &~~+\mu_2\left(\int_{-B}^{0}{\rm d}x~\rho_1(x)+\lambda-1+c \right)+\mu_1 \left(\int_{a}^{b}{\rm d}x~\rho_2(x)-c \right)\nonumber , \nonumber \end{align} \tag{ 141 }$$

where $\mu_1=A_1+A_2$ and $\mu_2=A_1$. Now minimising $\Sigma_c[\rho_I]$ with respect to $\rho_1$ and $\rho_2$, we get the following equations

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle &\frac{1}{2}x^2 - 2\alpha \int_{-B}^{0} {\rm d}y\rho_1(y)|x-y| -2 \alpha \lambda |x| -2 \alpha \int_{a}^{b} {\rm d}y~\rho_2(y)|x-y| + \mu_2=0, \nonumber \\ &~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~{\rm for}~-B\leqslant x <0, \label{eq-rho-122} \nonumber \end{align} \tag{ 142 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle &\frac{1}{2}x^2 - 2\alpha \int_{-B}^{0} {\rm d}y\rho_1(y)|x-y| -2 \alpha \lambda |x| -2 \alpha \int_{a}^{b} {\rm d}y~\rho_2(y)|x-y| + \mu_1=0, \nonumber \\ &~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~{\rm for}~a\leqslant x \leqslant b. \nonumber \end{align} \tag{ 143 }$$

Zoom In Zoom Out Reset image size

Taking derivative of the above two equations with respect to x on both sides and using $\newcommand{\f}{f} \frac{{\rm d}^2}{{\rm d}x^2}|x-y|=2 \delta(x-y)$, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \rho_1(x)=\frac{1}{4 \alpha},~~{\rm for}~-B\leqslant x <0, \label{rho-1} \nonumber \end{align} \tag{ 144 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \rho_2(x)=\frac{1}{4 \alpha},~~{\rm for}~a\leqslant x \leqslant b. \label{rho-2} \nonumber \end{align} \tag{ 145 }$$

We now insert the expression of $\rho_1(x)$ and $\rho_2(x)$ into (142), to find the following equation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle x\left [\frac{b-a}{2}+2\alpha \lambda -\frac{B}{2}\right] + \left[\mu_2-\frac{B^2}{4} - \frac{b^2-a^2}{4} \right] = 0 , \nonumber \end{align} \tag{ 146 }$$

which is valid for all x in the range $-B\leqslant x <0$. As a result we require that the coefficients of x and x⁰ in the above equation are zero. This implies the following two equations

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle 2 \alpha \lambda + \frac{(b-a)}{2} - \frac{B}{2}=0,~~\label{eq-B} \nonumber \end{align} \tag{ 147 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \mu_2 = \frac{B^2}{4}+\frac{b^2-a^2}{4} .~~~\label{mu_2} \nonumber \end{align} \tag{ 148 }$$

Similarly, inserting the expressions of $\rho_1(x)$ and $\rho_2(x)$ in (6.2) we find

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle 2 \alpha \lambda - \frac{(b+a)}{2} + \frac{B}{2}=0,~~\label{eq-b} \nonumber \end{align} \tag{ 149 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\f}{f} \displaystyle \mu_1 = \frac{B^2}{4}+\frac{b^2+a^2}{4} .~~~\label{mu_1} \nonumber \end{align} \tag{ 150 }$$

We have 6 unknowns ($a, ~b, ~B, ~\lambda, ~\mu_1, ~\mu_2$) to determine and till now we have 4 equations (147)–(150). We need two additional equations which are obtained from normalisations (see the last line of equation (6.2)): $\int_{-B}^0{\rm d}x~\rho_1(x)=1-c-\lambda$ and $\int_{a}^b{\rm d}x~\rho_2(x)=c$. This yields

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \frac{B}{4\alpha}=1-c-\lambda,~\label{eq-lam} \nonumber \end{align} \tag{ 151 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \frac{b-a}{4\alpha}=c.~~\label{eq-ab} \nonumber \end{align} \tag{ 152 }$$

Solving this system of six equations (147)–(152) for the six unknowns ($a, ~b, ~B, ~\lambda, ~\mu_1, ~\mu_2$), we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle &B=2\alpha, \nonumber \\ &b=2\alpha, \nonumber \\ &a= 2\alpha(1-2c), \nonumber \\ &\lambda=\frac{(1-2c)}{2}, \label{constants}\nonumber \\ &\mu_1=\alpha^2 (2+(1-2c)^2), \nonumber \\ &\mu_2=\alpha^2 (2-(1-2c)^2) . \nonumber \end{align} \tag{ 153 }$$

With these constants the equilibrium density in (140) is fully specified for $0\leqslant c \leqslant 1/2$ as (see equations (140), (144) and (145))

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \rho_{I}^*(x) = \left\{\begin{array}{@{}lcl@{}} \frac{1}{4\alpha} + \frac{(1-2c)}{2}~\delta(x),&{\rm for}&-2\alpha\leqslant x \leqslant 0 \nonumber \\ \frac{1}{4\alpha}&{\rm for}&2\alpha(1-2c)\leqslant x \leqslant 2\alpha . \end{array}\right. \label{den-eq} \nonumber \end{align} \tag{ 154 }$$

Finally inserting this expression of $\rho_{I}^*(x)$ in (141), we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \Sigma_c[\rho_{I}^*] = \frac{8 \alpha^2}{3}~(1/2-c)^3 -\frac{2\alpha^2}{3},~~{\rm for}~0\leqslant c \leqslant 1/2. \label{Sigma_c_lt_half} \nonumber \end{align} \tag{ 155 }$$

A similar computation for $1/2\leqslant c \leqslant 1$ yields

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \Sigma_c[\rho_{I}^*] = \frac{8 \alpha^2}{3}~(c-1/2)^3 -\frac{2\alpha^2}{3},~~{\rm for}~1/2\leqslant c \leqslant 1. \label{Sigma_c_gt_half} \nonumber \end{align} \tag{ 156 }$$

Combining both expressions we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \Sigma_c[\rho_{I}^*] = \frac{8 \alpha^2}{3}~|c-1/2|^3 -\frac{2\alpha^2}{3},~~{\rm for}~0\leqslant c \leqslant 1. \label{Sigma_c} \nonumber \end{align} \tag{ 157 }$$

Hence from equations (138) and (157), we have $\newcommand{\f}{f} \newcommand{\al}{\alpha} I(cN, N) \asymp {\rm e}^{-N^3[(8 \alpha^2/3)~|c-1/2|^3 -\frac{2\alpha^2}{3}]}$ to leading order for large N. Since $Z_N \approx I(N/2, N)$ and $N! \sim {\rm e}^{N \log N + O(N)}$, one finally obtains the large deviation form of the index distribution announced in equation (158)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle P_I(N_+=cN,N)\asymp \exp \left(-N^3\frac{8 \alpha^2}{3}~|c-1/2|^3 \right). \label{LDF} \nonumber \end{align} \tag{ 158 }$$

It is straightforward to check that this large deviation tail matches smoothly with the tails of the central region given in (134).

In this paper we have studied analytically the distribution of the position of the rightmost particle $x_{\max}$ of a 1d Coulomb gas confined in an external harmonic potential (the 1d jellium model) in the limit of large number of particles N. We have shown (see the appendix) that at very high temperature (when the inverse temperature $\beta = O(1)$) the distribution of $x_{\max}$, properly centred and scaled, is given by the Gumbel distribution. In the opposite limit of very low temperature (when $\newcommand{\g}{g} \beta = O(N^{1+\gamma})$ with $\newcommand{\g}{g} \gamma >0$), the distribution of $x_{\max}$ is given a given by a Gaussian centred at $\newcommand{\al}{\alpha} 2 \al$ with a width  ∼$\newcommand{\g}{g} N^{-1-\gamma/2}$ for large N. We have also shown that the most interesting case is when $\beta = O(N)$. In this case, we have shown that the cumulative distribution of $x_{\max}$, for large N, converges to a non-trivial limiting form $\newcommand{\al}{\alpha} F_\al(x)$ that is different from the Tracy–Widom distribution of the log-gas. We have shown that this function $\newcommand{\al}{\alpha} F_\al(x)$ is the solution of a non-local eigenvalue equation (17). We have also computed the rate functions associated with atypically large fluctuations around the mean (see equations (16), (19) and (20)) and found a third order phase transition between a pushed and a pulled phase, as in the log-gas.

In addition, we have studied the distribution of two other observables: (i) the gap $g=x_N-x_{N-1}$ between the two rightmost charges and (ii) the index N₊ which is the number of particles on the positive semi-axis. We have analytically computed the distribution of both quantities and found that their typical distributions can be expressed in terms of the same function $\newcommand{\al}{\alpha} F_\al(x)$ (see equations (95) and (100) for g and equations (132) and (133) for N₊). For both observables, the obtained limiting distributions are quite different from their counterpart in the log-gas. In both cases, we have computed the large deviations, to leading order for large N (see equations (103) and (158) for the gap and the index respectively).

Our work raises several interesting questions. For instance, how universal is the limiting distribution of $x_{\max}$ if one changes the confining potential or the pairwise repulsive interaction? It would be challenging to study $x_{\max}$ with a repulsive interaction of the form $|x_i-x_j|^{-k}$ (where $k \to 0$ corresponds to log-gas, while k  =  −1 corresponds to the 'jellium' model). Unlike the log-gas, the 1d jellium does not have a determinantal structure and computing its n-point correlations would be interesting.

The authors would like to acknowledge the support from the Indo-French Centre for the promotion of advanced research (IFCPAR) under Project No. 5604-2. AK would like to acknowledge the financial support from CNRS, France during his visit to LPTMS, Univ. Paris-Sud, Orsay where the paper has been finalised. This work was partially supported by ANR grant ANR-17-CE30-0027-01 RaMaTraF.

As mentioned in the text, the model can be studied at finite β, and the distributions of different observables can be computed for all β. In this appendix, we briefly describe the results for the distribution of $x_{\max}$ in the two limiting cases, namely very low and very high β.

We start with the very high temperature regime where $\beta = O(1)$. In this case, one can show that the interaction term plays a sub-dominant role and the particles essentially behave in an independent way, each of them sitting in a confining harmonic potential at thermal equilibrium at inverse temperature β. As a result the problem of finding the distribution of the position of the rightmost particle reduces effectively to computing the distribution of the maximum of N independent and identically distributed Gaussian random variables, each with zero mean. The limiting distribution for large N, properly centred and scaled, is thus given by the Gumbel distribution. More precisely, we find that for large N the typical fluctuations of $x_{\max}$ can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \displaystyle \label{scaling_gumbel} x_{\max} \approx 2\al - \sqrt{\frac{2 \ln N}{N}} + \frac{1}{\sqrt{2 N \, \ln N}} \, G \nonumber \end{align} \tag{ A.1 }$$

where G is a Gumbel random variable, i.e. ${\rm Prob.}(G \leqslant x) = {\rm e}^{-{\rm e}^{-x}}$.

In the opposite very low temperature limit, where $\newcommand{\g}{g} \beta = O(N^{1 +\gamma})$ with $\newcommand{\g}{g} \gamma >0$, the charges fluctuate independently around their equilibrium positions $\newcommand{\al}{\alpha} x_i^* = 2\al/N(2i-N-1)$ and the fluctuations are independent Gaussian around their respective means. Hence, in this case, computing the distribution of $x_{\max}$ again reduces to finding the maximum of a set of N independent random Gaussian variables, with the important difference that these random variables are no longer identical but they are Gaussians with mean at $x_i^*$ that depends explicitly on i. Consequently, the maximum is given precisely by x_N itself and hence $x_{\max}$ also has a Gaussian distribution. More precisely, we find that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\al}{\alpha} \newcommand{\f}{f} \newcommand{\g}{g} \displaystyle \label{scaling_Gauss} x_{\max} \approx 2 \al + \frac{1}{N^{1 + \frac{\gamma}{2}}} \, {{\mathcal N}}(0,1) , \nonumber \end{align} \tag{ A.2 }$$

where ${{\mathcal N}}(0, 1)$ is a zero mean, unit variance Gaussian random variable.