Exact extremal statistics in the classical $1d$ Coulomb gas

TL;DR

This paper analytically derives the probability distribution and large deviation functions for the maximum charge position in a 1D Coulomb gas, revealing a third-order phase transition and distinct fluctuations from the Tracy-Widom law.

Contribution

It provides the first exact analytical computation of the extremal statistics in the 1D Coulomb gas, including the distribution and phase transition analysis.

Findings

Typical fluctuations follow a non-trivial asymmetric scaling function.

The distribution differs from the Tracy-Widom law for Dyson's log-gas.

The system exhibits a third-order phase transition in the large deviation regime.

Abstract

We consider a one-dimensional classical Coulomb gas of $N$ like-charges in a harmonic potential -- also known as the one-dimensional one-component plasma (1dOCP). We compute analytically the probability distribution of the position $x_{\max}$ of the rightmost charge in the limit of large $N$. We show that the typical fluctuations of $x_{\max}$ around its mean are described by a non-trivial scaling function, with asymmetric tails. This distribution is different from the Tracy-Widom distribution of $x_{\max}$ for the Dyson's log-gas. We also compute the large deviation functions of $x_{\max}$ explicitly and show that the system exhibits a third-order phase transition, as in the log-gas. Our theoretical predictions are verified numerically.

Equations

$$\mathcal{P}(\{x_i\})=B_N{\mathrm{e}}^{-\frac{\beta }{2}\left(N{\sum }_{i=1}^N{x}_i^2-{\sum }_{i\ne j}\log (|x_i-x_j|)\right)},$$

$$Q(w,N)\approx {\mathcal{F}}_{\beta }\left(\sqrt{2}{N}^{2/3}(w-\sqrt{2})\right),$$

$$E(\{x_i\})=\frac{N^2}{2}\sum _{i=1}^N{x}_i^2-\alpha N\sum _{i\ne j}|x_i-x_j|,$$

$$\mathcal{P}(\{x_i\})=\frac{1}{Z_N}\exp [-E(\{x_i\})],$$

$$Q(w,N)=\text{Prob}.[x_{\max }\le w,N].$$

$$Q(w,N)\approx \{\begin{array}{ll}{\mathrm{e}}^{-N^3{\Phi }_{-}(w)+O(N^2)}& 0<2\alpha -w\sim O(1)\\ F_{\alpha }\left[N(w-2\alpha )+2\alpha \right]& |2\alpha -w|\sim O(1/N)\\ 1-e^{-N^2{\Phi }_{+}(w)+O(N)}& 0<w-2\alpha \sim O(1)\end{array}$$

$$\frac{d{F}_{\alpha }(x)}{dx}=A(\alpha ){\mathrm{e}}^{-x^2/2}{F}_{\alpha }(x+4\alpha ),$$

$$F_{\alpha }^{\prime }(x)\approx \{\begin{array}{ll}& \exp \left[-|x{|}^3/24\alpha +O(x^2)\right]\text{as}x\to -\infty \\ & \text{vspace}\ast -0.25cm\\ & \exp \left[-x^2/2+O(x)\right]\text{as}x\to \infty .\end{array}$$

$${\Phi }_{-}(w)=\{\begin{array}{ll}& \frac{(2\alpha -w{)}^3}{24\alpha },-2\alpha \le w\le 2\alpha \\ & \text{vspace}\ast -0.3cm\\ & \frac{w^2}{2}+\frac{2}{3}{\alpha }^2,w\le -2\alpha .\end{array}$$

$${\Phi }_{+}(w)$$

$$Q(w,N)$$

$$Z_N(w)$$

$$D_{\alpha }(x,N)=$$

$$Q(w,N)=\frac{D_{\alpha }(x,N)}{D_{\alpha }(\infty ,N)}\equiv F_{\alpha }(x,N).$$

$$\frac{d{F}_{\alpha }(x,N)}{dx}=\frac{D_{\alpha }(\infty ,N-1)}{D_{\alpha }(\infty ,N)}{\mathrm{e}}^{-\frac{x^2}{2}}{F}_{\alpha }(x+4\alpha ,N-1).$$

$${\rho }_w(x)=\frac{1}{4\alpha }+C\delta (x-w),-B\le x\le w,$$

$$C=1/2-w/(4\alpha ),B=2\alpha ,$$

$$P(x_{\max}=w,N)=\partial_{w}\,Q(w,N)\approx\left\{\begin{array}[]{rl}\mathrm{e}^{-N^{3}\,\Phi_{-}(w)}&\textrm{for}\,\,0<2\alpha-w\sim O(1)\;,\\ \\ N\,F_{\alpha}^{\prime}\left[N(w-2\alpha)+2\alpha\right]&\textrm{for}\,\,|2\alpha-w|\sim O(1/N)\\ \\ \mathrm{e}^{-N^{2}\,\Phi_{+}(w)}&\textrm{for}\,\,0<w-2\alpha\sim O(1)\;.\end{array}\right.$$

$$Q(w,N)$$

$$Z_N(w)$$

$$E(\{x_i\})$$

$${\rho }_w(x)=\frac{1}{N}\sum _{i=1}^N\delta (x-x_i).$$

$$E[\{x_i\}]=\mathcal{E}[{\rho }_w(x)]$$

$$\mathcal{E}[{\rho }_w(x)]=N^3\left(\frac{1}{2}{\int }_{-\infty }^{w}dx{x}^2{\rho }_w(x)-\alpha {\int }_{-\infty }^{w}dx{\int }_{-\infty }^{w}dy{\rho }_w(x){\rho }_w(y)|x-y|\right).$$

$$Z_N(w)=\int \mathcal{D}{\rho }_w\text{exp}\left(-\mathcal{E}[{\rho }_w(x)]\right)\delta \left({\int }_{-\infty }^{w}dx{\rho }_w(x)-1\right).$$

$$Z_N(w)$$

$$\mathcal{S}[{\rho }_w(x)]$$

$$Z_N(w)\approx \text{exp}\left(-N^{3}S[{\rho }_w^{\ast }(x)]\right),$$

$$\frac{1}{2}{x}^2-2\alpha {\int }_{-\infty }^{w}dy{\rho }_w^{\ast }(y)|x-y|+\mu =0.$$

$${\int }_{-B}^w{\rho }_w^{\ast }(x)dx=1.$$

$${\rho }_w^{\ast }(x)=\frac{1}{4\alpha }\mathrm{for}-2\alpha \le x\le 2\alpha ,w>2\alpha .$$

Full text

The Tracy-Widom (TW) distribution has emerged ubiquitously in diverse systems in the recent past buchanan ; quanta . It was originally discovered as the limiting distribution of the top eigenvalue $x_{\max}$ of an $N\times N$ Gaussian random matrix Tracy94 . Since then, it has appeared in various areas of physics Satya_lectures ; Krug2010 , mathematics BDJ1999 ; Joh2001 , and information theory Moustakas . For example, in physics it has appeared in stochastic growth models and related directed polymer in $1+1$ dimensional random media belonging to the Kardar-Parisi-Zhang (KPZ) universality class Joha2000 ; PS2000 ; MS2004 ; SS2010 ; CDR2010 ; Dot2010 ; ACQ2011 , non-intersecting Brownian motions FMS2011 , non-interacting fermions in a one-dimensional trapping potential Eisler2013 ; Dean15 ; Dean16 , disordered mesoscopic systems Beenaker and even in the Yang-Mills gauge theory in $2$-dimensions FMS2011 . It has also been measured experimentally in several systems including liquid crystals Takeuchi , coupled fiber lasers Davidson or disordered superconductors Lem2013 . The TW distribution describes the probability of typical fluctuations of $x_{\max}$ around its mean. In contrast, the atypical fluctuations of $x_{\max}$ to the left and right, far from its mean, are described respectively by the left and right large deviation tails. These tails have been computed explicitly BDG2001 ; Dean06 ; Dean08 ; Majumdar09 and shown to correspond to two different thermodynamic phases separated by a third order phase transition Nadal ; Majumdar14 . Similar third order phase transitions have also been found in a variety of other systems SMCF ; CP2013 ; Majumdar14 ; ATW2014 ; Ledou16 ; CMV2016 ; Baruch17 .

For Gaussian ensembles in random matrix theory (RMT), the joint probability distribution function (PDF) of the $N$ real eigenvalues $\{x_{1},\cdots,x_{N}\}$ is known explicitly Mehta91 ; Forrester10

[TABLE]

where $B_{N}$ is a normalisation constant and $\beta=1$, $2$ and $4$ depending on the symmetries of the matrices Mehta91 ; Forrester10 . This joint PDF can be interpreted as the equilibrium Gibbs distribution of a gas of $N$ charges with positions $x_{i}$’s that are confined on a line in the presence of an external harmonic potential and repelling each other via two dimensional logarithmic Coulomb interactions. This system is known as Dyson’s log-gas Dyson62 . In this picture the largest eigenvalue $x_{\max}=\max\{x_{1},\cdots,x_{N}\}$ corresponds to the position of the rightmost particle. The average density of eigenvalues $\rho_{N}(x)$ converges for large $N$ to the Wigner semi-circular law, $\rho_{\infty}(x)\sim\sqrt{2-x^{2}}/\pi$ which has a finite support $[-\sqrt{2},+\sqrt{2}]$. Hence the average $\langle x_{\max}\rangle\sim\sqrt{2}$ for large $N$. The typical fluctuations of $x_{\max}$ around its mean $\sqrt{2}$ are of ${O}(N^{-2/3})$. On this scale the cumulative distribution $Q(w,N)={\rm Prob}(x_{\max}\leq w,N)$, takes the scaling form

[TABLE]

where ${\cal F}_{\beta}(x)$ is the Tracy-Widom distribution. This scaling function can be written in terms of the solution of a Painlevé II equation Tracy94 and has non-Gaussian tails. Interestingly, even though the TW distribution was derived originally for a harmonic confining potential, it has been shown to be universal with respect to the shape of the confining potential, as long as the average density vanishes at the upper edge as a square root (as in the case of the Wigner semi-circular law). A natural question then arises whether the TW distribution for $x_{\max}$ is robust when one changes, instead of the confining potential, the nature of the repulsive interaction between the charges.

A natural candidate model to address this question is the system of one-dimensional charges in a harmonic potential but interacting pairwise via the true $1d$ Coulomb potential. The energy of this system is given by

[TABLE]

where $\alpha\geq 0$ denotes the strength of the Coulomb repulsion. The choice of the powers of $N$ in the coupling constants is such that for large $N$ (i) both terms in the energy are of same order and (ii) the charges are confined in a region whose span is $O(1)$. Indeed this is also the 1dOCP or “jellium” model, where $N$ charges of the same sign interact in the presence of a uniform background of opposite charges, assuring charge neutrality. This model is a paradigm for $1d$ charged plasma Choquard1981 as several observables can be calculated analytically Lenard61 ; Prager62 ; Baxter63 ; Dean2010 ; TT2015 .

In this Letter, we show that the statistics of the position of the rightmost charge $x_{\max}$ can also be computed exactly in this 1dOCP model. Our main result is to show that the limiting distribution of the typical fluctuations of $x_{\max}$ in this model is indeed different from the TW distribution, found in the Dyson’s log-gas. Moreover, by computing the left and the right large deviation functions explicitly, we find that the third-order phase transition between a pushed gas (left large deviation) and a pulled gas (right large deviation) is still present in the 1dOCP model, as in the case of the log-gas. Interestingly, a similar third-order phase transition between the pushed and the pulled phase was recently found Cunden2017 by analysing large deviation functions associated with the position of the farthest charge in a $d$-dimensional jellium model, though the limiting distribution of the position of the farthest charge is still open for this $d$-dimensional problem. In $d=1$ this corresponds to the distribution of the maximum of $|x_{i}|$’s of the charges.

We start with the joint PDF of the positions $x_{i}\in(-\infty,\infty)$ in the 1dOCP, given by the Boltzmann weight

[TABLE]

where $Z_{N}$ is the partition function and the energy $E(\{x_{i}\})$ is given in (3). In the large $N$ limit, the average density can be obtained by minimizing the energy $E(\{x_{i}\})$. It is easy to show that the minimum energy configuration is given by $x_{j}^{*}=\frac{2\alpha}{N}(2j-N-1)$ ($j=1,\cdots,N$). Hence the particles are equi-spaced and the rightmost (leftmost) particle is at $x_{N}^{*}=2\alpha(1-1/N)$ (respectively at $x_{1}^{*}=-2\alpha(1-1/N)$). This implies that in the $N\to\infty$ limit, the average density profile $\rho_{\infty}(x)$ is flat: $\rho_{\infty}(x)=\frac{1}{4\alpha}$ for $-2\alpha\leq x\leq 2\alpha$ (see Fig. 1), in contrast to the Wigner semi-circle in the log-gas. Our focus here is on the large $N$ behavior of the cumulative distribution of the rightmost particle,

[TABLE]

To anticipate the scaling behavior of $Q(w,N)$, we first make the following observations. It follows from the above analysis of the average density that the mean position of the rightmost particle is at $\langle x_{\max}\rangle=x_{N}^{*}=2\alpha(1-1/N)$. Given that the average density is uniform, the typical separation between two adjacent particles is of $O(1/N)$ everywhere. Hence, the scale of typical fluctuations of $x_{\max}$ around its average is also of $O(1/N)$. This suggests that the cumulative probability distribution $Q(w,N)$, in the region of typical fluctuations where $|w-2\alpha|\sim O(1/N)$, should have the scaling form for large $N$, $Q(w,N)\approx F_{\alpha}(N(w-x_{N}^{*}))=F_{\alpha}(N(w-2\alpha)+2\alpha)$ where $F_{\alpha}(x)$ is a nontrivial scaling function (the analogue of the TW distribution for the log-gas). In this paper, we compute this scaling function $F_{\alpha}(x)$. In addition, for atypical fluctuations where $x_{\max}-\langle x_{\max}\rangle\sim O(1)$, both to the left and to the right of the mean, $Q(w,N)$ has large deviation tails that are also computed explicitly. More precisely, we find

[TABLE]

where $\Phi_{-}(w)$ and $\Phi_{+}(w)$ are the left and right rate functions. We show that the scaling function $F_{\alpha}(x)$ in the central regime satisfies a nonlocal eigenvalue equation

[TABLE]

with the boundary conditions: $F_{\alpha}(-\infty)=0$ and $F_{\alpha}(\infty)=1$. These boundary conditions, along with the fact that $F_{\alpha}(x)\geq 0$ for all $x$, uniquely fixes the eigenvalue $A(\alpha)$. Clearly, the scaling function $F_{\alpha}(x)$ is different from the TW distribution. While it is hard to compute $A(\alpha)$ explicitly for all $\alpha\geq 0$ (for a numerical plot of $A(\alpha)$, see Fig. 2), we can determine its small and large $\alpha$ behaviors: $A(\alpha)\to 1/(4\,e\,\alpha)$ as $\alpha\to 0$ and $A(\alpha)\to 1/\sqrt{2\pi}$ as $\alpha\to\infty$. From Eq. (7), we can derive the leading asymptotic tails of the PDF $F_{\alpha}^{\prime}(x)$ for all $\alpha$

[TABLE]

We note that the leading $x\to-\infty$ behavior of $F_{\alpha}^{\prime}(x)$ is identical to that of the TW distribution ${\cal F}^{\prime}_{\beta=1/\alpha}(x)$, while the right tail of $F_{\alpha}^{\prime}(x)$ decays faster than the TW right tail Dumaz2013 . This is indeed our main result. In addition, we also compute exactly the large deviations rate functions. For the left tail we find

[TABLE]

For the right tail, we find

[TABLE]

It is easy to check, using the asymptotic behavior of $F_{\alpha}^{\prime}(x)$ for large $|x|$ in Eq. (8), that the central part of the distribution of $x_{\max}$ matches smoothly with the two large deviation regimes flanking this central part. Indeed, as discussed later, the vanishing of $\Phi_{-}(w)$ when $w\to 2\alpha$ as a cubic power is responsible for a third order phase transition at the critical point $w=2\alpha$, in very much the same way as in the log-gas Majumdar14 .

We start from the joint PDF of $\{x_{i}\}$’s in Eq. (4). We note that $Q(w,N)=\text{Prob}.[x_{\max}\leq w]=\text{Prob.}(x_{1}\leq w,\cdots,x_{N}\leq w)$. Hence it can be expressed as the ratio of two partition functions

[TABLE]

with $E(\{x_{i}\})$ given in (3) and we have suppressed the $\alpha$-dependence in $Z_{N}(w)$ for simplicity. Note that $Z_{N}(w)$ can be interpreted as the partition function of the 1dOCP in presence of a hard wall at $w$. Below, we analyse $Q(w,N)$ in the central regime first, followed by the two large deviation tails.

Central regime: Noting that the energy function $E(\{x_{i}\})$ in (3) is symmetric under permutations over the positions $(x_{1},x_{2},\cdots,x_{N})$, we write

[TABLE]

where $\Theta(x)$ is the Heaviside theta function. For an ordered configuration $(x_{1}<x_{2}<\cdots<x_{N})$, one can eliminate the absolute values $|x_{i}-x_{j}|$ ’s and rewrite the energy function $E(\{x_{i}\})$ in (3) as $E(\{x_{i}\})=\frac{N^{2}}{2}\sum_{i=1}^{N}\left(x_{i}-\frac{2\alpha}{N}(2i-N-1)\right)^{2}+C_{N}(\alpha)$, where $C_{N}(\alpha)=2\alpha^{2}\sum_{i=1}^{N}\left(2i-N-1\right)^{2}$ is just a constant. This trick of eliminating the absolute values via ordering has been used before for 1dOCP in numerous contexts Lenard61 ; Baxter63 ; Dean2010 ; TT2015 . Performing a change of variables $\epsilon_{k}=\left(Nx_{k}-2\alpha(2k-N-1)\right)$ for all $k=1,2,\cdots,N$ in (13), we can rewrite $Z_{N}(w)=N!~{}D_{\alpha}\left(N\left(w-\frac{2\alpha}{N}(N-1)\right),N\right)$ where

[TABLE]

Therefore setting $x=N\left(w-\frac{2\alpha}{N}(N-1)\right)$, $Q(w,N)$ in Eq. (20) can be written as

[TABLE]

Taking derivative with respect to $x$ in Eq. (15), and using Eq. (14), we obtain

[TABLE]

To estimate the ratio ${D_{\alpha}(\infty,N-1)}/{D_{\alpha}(\infty,N)}$ for large $N$, we note from Eq. (14) that $D_{\alpha}(\infty,N)$ can be interpreted as the partition function of an auxiliary gas of particles with positions $\xi_{1},\xi_{2},\cdots,\xi_{N}$ confined in an external harmonic potential and with the one-sided constraint $\xi_{k-1}<\xi_{k}+4\alpha$ for all $k=2,3,\cdots N$. Indeed this constraint provides a short-range interaction between the particles. Thus our original problem of the 1dOCP which has long-range interaction is mapped onto a problem of short-ranged interacting particles. For such a short-ranged system, it is natural to expect that the free energy is extensive in $N$. Thus one would expect that, for large $N$, the partition function must scale as $D_{\alpha}(\infty,N)\sim[A(\alpha)]^{-N}$ where $\ln A(\alpha)$ is the free energy per particle. Thus the ratio ${D_{\alpha}(\infty,N-1)}/{D_{\alpha}(\infty,N)}\to A(\alpha)$ as $N\to\infty$. This suggests that $F_{\alpha}(x,N)$ should converge to a limiting form $F_{\alpha}(x)$ for large $N$, which then satisfies the nonlocal eigenvalue equation (7). Thus the eigenvalue $A(\alpha)$ has a physical interpretation: $\ln A(\alpha)$ is the free energy per particle of a short-ranged interacting gas. However, computing analytically $A(\alpha)$ for all $\alpha$ seems hard. Interestingly, Baxter Baxter63 encountered a similar nonlocal eigenvalue equation while computing the partition function of the 1dOCP in a finite box $[-L,+L]$ and analysed the eigenvalue $A(\alpha)$ in the two limits $\alpha\to 0$ and $\alpha\to\infty$. Translating his results to our problem, following a simple rescaling of the parameters, we obtain the asymptotic results for $A(\alpha)$ announced before.

It is straightforward to derive the asymptotic tails of $F_{\alpha}^{\prime}(x)$ in Eq. (8). We consider first the limit $x\to\infty$ where $F_{\alpha}(x+4\alpha)\to 1$ on the right hand side (rhs) of Eq. (8). Hence, to leading order, $F_{\alpha}^{\prime}(x)\approx A(\alpha)\,\mathrm{e}^{-x^{2}/2}$, providing the Gaussian right tail in Eq. (8). To compute the left tail, we make a simple ansatz that $F_{\alpha}(x)\approx\mathrm{e}^{-a_{0}\,|x|^{\delta}}$ as $x\to-\infty$, where $a_{0}$ and $\delta$ are to be determined. Substituting this ansatz in the rhs of Eq. (8) yields $\approx A(\alpha)\,\mathrm{e}^{-x^{2}/2-a_{0}\,(|x|-4\alpha)^{\delta}}$. For large $|x|$, $(|x|-4\alpha)^{\delta}\sim|x|^{\delta}(1-4\,\alpha\,\delta/|x|)$ to leading orders. Hence the rhs behaves as $A(\alpha)\mathrm{e}^{-a_{0}\,|x|^{\delta}-x^{2}/2+4\,\alpha\,\delta\,|x|^{\delta-1}}$. The left hand side (lhs) of Eq. (8) behaves as $\approx\mathrm{e}^{-a_{0}\,|x|^{\delta}}$ to leading order. Comparing both sides, we see that the term $x^{2}/2$ and $|x|^{\delta-1}$ on the rhs must cancel each other implying $\delta=3$ and $a_{0}=1/(8\alpha\,\delta)=1/(24\,\alpha)$. This provides the leading left tail of $F_{\alpha}^{\prime}(x)$ in Eq. (8).

For general $\alpha>0$, determining explicitly the eigenvalue $A(\alpha)$ and the full scaling function $F_{\alpha}(x)$ is difficult. However they can be obtained by solving (7) numerically by tuning the value of $A(\alpha)$ using the standard shooting method. This gives $F_{\alpha}(x)$ and $A(\alpha)$ simultaneously. In Fig. 2 (left panel), we plot $A(\alpha)$ vs. $\alpha$ and compare with its predicted asymptotics. In Fig. 2 (right panel), we compare $F_{\alpha}^{\prime}(x)$ evaluated numerically from this shooting method, with the one obtained from direct Monte-Carlo simulation of the 1dOCP. The agreement is excellent.

Left large deviation function. We consider $Q(w,N)$ in Eqs. (20) and (12) with $0<2\alpha-w\sim O(1)$. Since $w$ represents the position of the hard wall, $w<2\alpha$ corresponds to “pushing” the charges to the left of the right edge at $2\alpha$. This disturbs the originally flat density and leads to a collective reorganisation of all the charges, as in the case of the log-gas Dean06 ; Dean08 . We get instead a new equilibrium density that minimizes the energy, i.e., a new saddle point of the integral in Eq. (12). It is well known that in the jellium model, the bulk density is insensitive to the location of a wall ASZ2014 . This implies that in the bulk, the density is still given by the original equilibrium value $1/(4\alpha)$, for $-2\alpha<w\leq 2\alpha$. Hence, when the wall moves to the left of $2\alpha$, all the charges that get pushed by the wall must get absorbed at the wall. This observation leads us to look for a saddle point density of the form

[TABLE]

where the constant bulk density is supported over the interval $[-B,w]$. We then minimize the energy with respect to the two variational parameters $B$ and $C$. Skipping details (see supp_mat ), we find that

[TABLE]

as long as $-2\alpha\leq w\leq 2\alpha$. When $w$ hits $-2\alpha$ from the right, all the charges get absorbed at the wall $w$ and the saddle point density is just $\rho_{w}(x)=\delta(x-w)$, for all $w\leq-2\alpha$. Substituting $\rho_{w}(x)$ in the energy supp_mat , we find the results for $\Phi_{-}(w)$ given in Eq. (9).

Right large deviations. For fluctuations $(x_{\max}-2\alpha)\sim{O}(1)$ to the right of the edge $2\alpha$, we consider the PDF $\partial_{w}Q(w,N)$ in Eqs. (20) and (12) with $w>2\alpha$. It turns out that the configuration that dominates this integral is one where the rightmost charge is at $w>2\alpha$, while the rest of the $N-1$ charges stay in the equilibrium configuration with a flat profile over the interval $[-2\alpha,+2\alpha]$. This is analogous to the “pulled” phase in the log-gas Majumdar09 . Thus, for large $N$, the PDF can be approximated as $\partial_{w}Q(w,N)\approx\mathrm{e}^{-\Delta E_{\rm pulled}}$, where $\Delta E_{\rm pulled}$ is the energy cost of pulling the rightmost particle from the “sea” of $N-1$ particles in the equilibrium flat configuration. This energy cost can be estimated from Eq. (3): a first contribution from the change in the external potential energy of the rightmost charge and a second due to the interaction of the rightmost charge with the $(N-1)$ other particles. One gets (see supp_mat for details): $\Delta\,E_{\rm pulled}\approx N^{2}\left(\frac{w^{2}}{2}-\frac{1}{2}\int_{-2\alpha}^{2\alpha}(w-x)\,dx\right)$ up to a constant. This gives $\partial_{w}Q_{N}(w)\approx\mathrm{e}^{-N^{2}\Phi_{+}(w)}$ where $\Phi_{+}(w)$ is given in (10).

Since $Q(w,N)$ is the ratio of two partition functions (12), $-\ln Q(w,N)$ is the free energy difference between the pushed (left) and the pulled (right) phase. From Eq. (6), this free energy $\propto\Phi_{-}(w)$ has a singular behavior at the critical point $w=2\alpha$. Indeed it vanishes as a cubic power as $w\to 2\alpha$ from the left [see Eq. (9)], leading to a discontinuity of the third derivative of $\Phi_{-}(w)$ at $w=2\alpha$. This third order phase transition at $w=2\alpha$ is similar to the one in the log-gas. Unlike in the log-gas, there is an additional third-order phase transition in this 1dOCP when $w\to-2\alpha$ [see Eq. (9)]. However, this transition is not of the “pushed-pulled” type like the one at $w=2\alpha$, but rather a condensation-type transition as all charges accumulate at the wall for $w\leq-2\alpha$.

Conclusion: In this Letter 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 (1dOCP), in the limit of large number of particles $N$. We have obtained the limiting large $N$ distribution describing the typical fluctuations of $x_{\max}$ around its mean, and shown that it is different from the Tracy-Widom distribution of the log-gas. We also computed the rate functions associated with atypically large fluctuations around the mean and found a third order phase transition between a pushed and a pulled phase, as in the log-gas. 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 1dOCP). Unlike the log-gas, the 1dOCP does not have a determinantal structure and computing its $n$-point correlations would be interesting.

Acknowledgements. We thank M. Krishnapur, D. Mukamel, E. Trizac and P. Vivo for discussions. We acknowledge support from the Indo- French Centre for the Promotion of Advanced Research (IFCPAR) under Project 5604-2.

I Left large deviation function

To compute the left large deviation function, we start with the cumulative probability $Q(w,N)=\textrm{Prob.}(x_{\max}\leq w,N)$ which is given by

[TABLE]

and $\alpha>0$. In the large $N$ limit, the leading behavior of the partition function $Z_{N}(w)$ of the Coulomb gas with a wall at $w$ can be computed as follows. We first introduce a macroscopic charge density, normalized to unity with its support on $(-\infty,w]$

[TABLE]

The energy of the gas of particles in (22), to leading order for large $N$, can be written in terms of the macroscopic density $\rho_{w}(x)$ as

[TABLE]

where

[TABLE]

The $N$-fold integration in the partition function $Z_{N}(w)$ in Eq. (21) can be carried out in two steps. In the first step, one fixes the macroscopic density $\rho_{w}(x)$ and one sums over all the microscopic configurations of $x_{i}$’s consistent with this density $\rho_{w}(x)$. In the second step, one sums over all possible macroscopic densities $\rho_{w}(x)$ that are positive and normalized to unity $\int_{-\infty}^{w}\rho_{w}(x)\,dx=1$. The first step gives rise to an entropy term that scales, for large $N$, as $\mathrm{e}^{O(N)}$. Since this is much smaller than the energy term, which scales like $\mathrm{e}^{O(N^{3})}$ from Eq. (25), we henceforth neglect the entropy term. Neglecting this entropy term, the partition function can then be expressed as a functional integral

[TABLE]

Replacing the delta function by its integral representation, we have

[TABLE]

For large $N$, we can now evaluate this functional integral by saddle point method. This gives

[TABLE]

where $\rho^{*}_{w}(x)$ is the saddle point density that minimizes the action $S[\rho_{w}(x)]$ in Eq. (28). The equation for $\rho^{*}_{w}(x)$ is obtained from $\left(\frac{\delta S[\rho_{w}]}{\delta\rho_{w}(x)}\right)_{\rho_{w}=\rho^{*}_{w}}=0$ as

[TABLE]

Let us now specify the limits of the integration in the above equation. This equation holds for $x$ belonging to the support of $\rho_{w}^{*}(x)$. Clearly the support can not be on $(-\infty,w]$. This can be seen easily from the observation that when $x\to-\infty$, the first term in Eq. (30) scales as $x^{2}$, while the second term scales as $|x|$ – hence they can not compensate each other. Therefore the support must be over a finite region $[-B,w]$ where $B$ is to be determined from the normalization condition

[TABLE]

For $x\in[-B,w]$, it is easy to see by differentiating twice the saddle point equation (30) and using the identity $\frac{d^{2}}{dx^{2}}|x-y|=2\delta(x-y)$, that $\rho^{*}_{w}(x)=1/(4\alpha)$. Clearly, if $w>2\alpha$, the saddle point density is given by

[TABLE]

Thus for $w>2\alpha$ 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 $w<2\alpha$, the wall tries to push the charges to the left of $2\alpha$ (see the left panel of Fig. 3). We have seen from above that the bulk density does not change from its equilibrium value $\rho^{*}_{w}(x)=1/(4\alpha)$ to the left of the wall at $w$. Normalization to unity of the charge density means that the extra charge that the wall displaces must be absorbed at the wall, since the bulk is not affected. This leads, for $w<2\alpha$, to a new saddle point density of the form

[TABLE]

where $C$ represents the density of the charges displaced and absorbed at the wall. The normalization condition $\int_{-B}^{w}\rho^{*}_{w}(x)\,dx=1$ relates the two parameters $B$ and $C$ via

[TABLE]

Substituting the saddle point density $\rho^{*}_{w}(x)$ in (30) yields

[TABLE]

Performing the integral over $y$ explicitly, we get

[TABLE]

Since Eq. (36) is valid for arbitrary $x\in[-B,w]$, we get two additional equations

[TABLE]

We therefore have three equations (34), (37) and (38) for three unknowns $\mu,B$ and $C$. Solving, we get

[TABLE]

Note that the condition $C\leq 1$ (needed for the normalization) indicates that the above analysis is only valid for $w>-2\alpha$. When $w\to-2\alpha$, $C\to 1$: this means that all the charges are absorbed at the wall and there is no bulk charge left. Thus for $w<-2\alpha$, 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$

[TABLE]

Let us first consider $w>2\alpha$. In this case $\rho_{w}^{*}(x)=1/(4\alpha)$ for $x\in[-2\alpha,+2\alpha]$. Substituting this density in Eq. (28) we get the saddle point action

[TABLE]

Therefore from Eq. (29) the partition function $Z_{N}(w)$ for large $N$ and for $w>2\alpha$ behaves as

[TABLE]

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

[TABLE]

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

[TABLE]

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 know consider the region where $-2\alpha\leq w\leq 2\alpha$. Substituting the saddle point density $\rho_{w}^{*}(x)$ from the second line of Eq. (42) in Eq. (28) we get

[TABLE]

Substituting this result in Eq. (29) and using the expression for the denominator in Eq. (45) we get

[TABLE]

Finally, we consider the region where $w\leq-2\alpha$. Substituting the saddle point density $\rho_{w}^{*}(x)$ from the third line of Eq. (42) in Eq. (28) we get

[TABLE]

Substituting this result in Eq. (29) and using the expression for the denominator in Eq. (45) we get

[TABLE]

Summarizing we obtain the result in Eq. (9) of the main text

[TABLE]

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

Third order phase transition at $w=2\alpha$. The cumulative distribution $Q(w,N)$ in Eq. (20) 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 Eq. (46), we see that to leading order for large $N$, $-\ln Q(w,N)\approx 0$ for $w>2\alpha$. In contrast, for $w<2\alpha$, using Eq.  (48) and (50), we see that $-\ln Q(w,N)\approx N^{3}\Phi_{-}(w)$ where $\Phi_{-}(w)$ is given in Eq. (51). Hence, we get (see Fig. 4)

[TABLE]

Thus $\Phi_{-}(w)$ is just the free energy cost in pushing the wall $w$ to the left of the right edge $2\alpha$ (see Fig. 3). From the expression of $\Phi_{-}(w)$ in the first line of Eq. (51), it follows that $\Phi_{-}(w)$ vanishes as the third power $\Phi_{-}(w)\propto(2\alpha-w)^{3}$ as $w\to 2\alpha$ from the left. Thus the third derivative of the free energy difference vanishes at the critical point $2\alpha$, making this a third order phase transition. Indeed, one can also look at the pressure on the wall, which is simply ${\cal P}=-N^{3}\,\Phi^{\prime}_{-}(w)$ (derivative of the free energy with respect to the wall position). Clearly, the pressure ${\cal P}$ is zero for $w>2\alpha$ (the charges do not touch the wall) and is non zero for $w<2\alpha$. The pressure ${\cal P}$ vanishes as $N^{3}\,(2\alpha-w)^{2}$ as $w\to 2\alpha$ from below.

II Right large deviation function

Here we focus, for large $N$, on the distribution $Q(w,N)$ in the region $0<w-2\alpha\sim O(1)$. From the analysis in the previous section, we have seen that in this regime, to leading order for large $N$, $Q(w,N)\approx 1$ [see Eq. (46)]. To compute the subleading corrections to this leading order term $1$, it is convenient to consider instead the PDF of $x_{\max}$, given by the derivative of Eq. (20)

[TABLE]

where we have set $x_{N}=w$ in Eqs. (21) and (22). This can be re-written as

[TABLE]

where $\langle\cdots\rangle_{N-1}$ denotes the average over the $N-1$ charges. We can then analyse this average for large $N$, for $w>2\alpha$, following Ref. sm-Majumdar09 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 $w>2\alpha$, while the rest of $N-1$ charges should be in their equilibrium flat configuration, i.e., with a density $\rho^{*}_{w}(x)=1/(4\alpha)$ for $x\in[-2\alpha,2\alpha]$ (see the right panel of Fig. 3). Furthermore, for large $N$, to leading order, we can approximate

[TABLE]

where $C_{0}$ is a constant, i.e. independent of $w$. Using $\rho^{*}_{w}(x)=1/(4\alpha)$ for $x\in[-2\alpha,2\alpha]$ and performing the integral in Eq. (55), we obtain

[TABLE]

where

[TABLE]

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