Harmonically confined particles with long-range repulsive interactions

TL;DR

This paper analyzes a classical particle system confined by a harmonic trap with long-range repulsive interactions, deriving exact density profiles for large N and revealing complex dependence on the interaction parameter k.

Contribution

It provides an exact computation of the average density profile for a broad class of long-range interacting particles under harmonic confinement, generalizing known models.

Findings

Density profile is independent of temperature at low temperatures.

Distinct behaviors of the density profile for different ranges of k.

Exact solutions for large N across all k > -2.

Abstract

We study an interacting system of $N$ classical particles on a line at thermal equilibrium. The particles are confined by a harmonic trap and repelling each other via pairwise interaction potential that behaves as a power law $\propto \sum_{\substack{i\neq j}}^N|x_i-x_j|^{-k}$ (with $k>-2$) of their mutual distance. This is a generalization of the well known cases of the one component plasma ($k=-1$), Dyson's log-gas ($k\to 0^+$), and the Calogero-Moser model ($k=2$). Due to the competition between harmonic confinement and pairwise repulsion, the particles spread over a finite region of space for all $k>-2$. We compute exactly the average density profile for large $N$ for all $k>-2$ and show that while it is independent of temperature for sufficiently low temperature, it has a rich and nontrivial dependence on $k$ with distinct behavior for $-2<k<1$, $k>1$ and $k=1$.

Equations

$$E(\{x_i\})=\frac{1}{2}\sum _{i=1}^N{x}_i^2+\frac{J\text{sgn}(k)}{2}\sum _{i\ne j}\frac{1}{|x_i-x_j{|}^k}$$

$$E[\{x_i\}]=\frac{1}{2}\sum _{i=1}^N{x}_i^2-\frac{1}{2}\sum _{i\ne j}\ln |x_i-x_j|.$$

$$⟨{\rho }_N(x)⟩\approx \frac{1}{\sqrt{N}}{\tilde{\rho }}_{\mathrm{sc}}\left(\frac{x}{\sqrt{N}}\right),{\tilde{\rho }}_{\mathrm{sc}}(y)=\frac{1}{\pi }\sqrt{2-y^2}.$$

$${\alpha }_k=\{\begin{array}{ll}& \frac{1}{k+2},-2<k<1\\ & \frac{k}{k+2},k>1.\end{array}$$

$$F_k(z)=\frac{1}{B({\gamma }_k+1,{\gamma }_k+1)}{\left(\frac{1}{4}-z^2\right)}^{{\gamma }_k}.$$

$${\gamma }_k=\{\begin{array}{ll}& \frac{k+1}{2},-2<k<1\\ & \frac{1}{k},k>1.\end{array}$$

$$Z_N(\beta )\sim \int d\mu \int \mathcal{D}[{\tilde{\rho }}_k]\exp \left(-\beta N^{\frac{4+k}{2+k}}\Sigma [{\tilde{\rho }}_k(y)]\right)$$

$$\Sigma [{\tilde{\rho }}_k(y)]$$

$$\frac{y^2}{2}+J\text{sgn}(k)\int d{y}^{\prime }\frac{{\tilde{\rho }}_k(y^{\prime })}{|y^{\prime }-y{|}^k}=\mu .$$

$$PV{\int }_{-{\ell }_k/2}^{{\ell }_k/2}\frac{\text{sgn}(y^{\prime }-y)}{|y-y^{\prime }{|}^{k+1}}{\tilde{\rho }}_k(y^{\prime })d{y}^{\prime }=-\frac{y}{J|k|},k\ne 0$$

$$\mathcal{E}\approx \frac{N}{2}\int dx{x}^2{\rho }_N(x)+J\zeta (k)N^{k+1}\int dx{\left[{\rho }_N(x)\right]}^{k+1},$$

$$Z_N(\beta )\sim \int d\mu \int \mathcal{D}[{\tilde{\rho }}_k]\exp \left(-\beta N^{\frac{3k+2}{k+2}}\Sigma [{\tilde{\rho }}_k(y)]\right)$$

$$\Sigma [{\tilde{\rho }}_k(y)]$$

$$\frac{1}{2}{y}^2+J\zeta (k)(k+1){\left[{\tilde{\rho }}_k(y)\right]}^k=\mu .$$

$${\tilde{\rho }}_1(y)=\frac{1}{4J}(2\mu -y^2),\mu =\frac{1}{2}(3J{)}^{2/3}.$$

Full text

Introduction: A gas of $N$ classical particles, confined by a harmonic potential on a line and interacting with each other via pairwise repulsion, constitutes one of the simplest interacting particle systems that have been well studied in the past. It has seen a recent revival in the wake of the physics of cold atoms. When the pairwise repulsive interaction decays as a power-law of the distance between the particles, the energy of the so called Riesz gas Riesz is given by

[TABLE]

where $J>0$ and $\{x_{i}\}$ ($i=1,2,\ldots,N$) denote the positions of the particles on the line. The index $k>-2$ characterizes the strength of the pairwise interaction and ${\rm{sgn}}(k)$ in the prefactor ensures a repulsive interaction. For $k<-2$, the quadratic potential is not strong enough to counter the strong repulsion and confine the particles. Consequently the particles fly off to $\pm\infty$ and thus the case $k<-2$ is not physically interesting. Given the energy in (1), the joint probability distribution function (PDF) of the particles’ positions is given by the Boltzmann weight, $P(x_{1},\cdots,x_{N})=\frac{1}{Z_{N}(\beta)}e^{-\beta E[\{x_{i}\}]}$, where $\beta$ is the inverse temperature ($\rm k_{B}=1$) and $Z_{N}(\beta)=\int\prod_{i=1}^{N}{dx_{i}}\,e^{-\beta E[\{x_{i}\}]}$ is the normalizing partition function. The harmonic potential tries to confine the particles near the center of the trap, while the repulsive interaction tries to push them apart. As a result of the competition between the two terms, it turns out that the particles get confined to a finite region of space for large $N$, with a space-dependent average macroscopic density, $\langle\rho_{N}(x)\rangle=N^{-1}\,\sum_{i=1}^{N}\langle\delta(x-x_{i})\rangle$ (normalized to unity), where $\langle\ldots\rangle$ denotes an average with respect to the Boltzmann weight. A basic natural question is: what is the configuration of $x_{i}$’s that minimizes the energy in (1) for large $N$ and what is the density profile in the ground state ? This is a classic and important optimization problem both in physics (see below) and in mathematics (see e.g. Refs. Land72, ; CGZ14, ; LS17, ) whose solution, for generic $k>-2$, is hitherto unknown. A related question is: how does the average density profile depend on the inverse temperature $\beta$ ?

This problem is of great general interest as there are varied physical systems that correspond to special values of $k$. We start with $k=-1$ where the interaction is linearly repulsive with distance. This is the well known one dimensional one-component plasma (1$d$OCP) Forrester_book , consisting of oppositely charged particles with pairwise Coulomb interaction (linear in $1d$) and overall charge neutrality Lenard_61 ; Prager_62 ; Baxter_63 ; Dhar et al. (2017); dhar2018 . Integrating out the positions of the negative charges gives rise to an effective quadratic confinement for the positive charges and the effective energy of the $N$ positive charges with coordinates $\{x_{i}\}$ is precisely given by (1) with $k=-1$. In this case, the energy can be easily minimized by ordering the positions of the particles leading to an equispaced configuration Lenard_61 ; Prager_62 ; Baxter_63 ; Dhar et al. (2017); dhar2018 . Moreover, for large $N$, the average density profile $\langle\rho_{N}(x)\rangle$ turns out to be independent of $\beta$ and approaches a scaling form $\langle\rho_{N}(x)\rangle\to(1/N)\tilde{\rho}_{\rm OCP}(x/N)$, where the scaled density $\tilde{\rho}_{\rm OCP}(y)=1/(2J)$ is uniform over the interval $[-J,+J]$ and vanishes outside CGZ14 ; Lenard_61 ; Prager_62 ; Baxter_63 ; Dhar et al. (2017); dhar2018 ; Cunden_2017 .

The second and perhaps the most well studied example corresponds to the limit $k\to 0^{+}$, where we replace ${\rm sgn}(k)$ in Eq. (1) by $+1$, use $|x_{i}-x_{j}|^{(-k)}\approx 1-k\log|x_{i}-x_{j}|$ and set $J=1/k$. The energy in (1) then reduces, up to an overall additive constant, to

[TABLE]

This is the celebrated log-gas of Dyson Dyson_1962 . For the special values of $\beta=1$, $2$ and $4$, the Boltzmann weight of the log-gas $P(\{x_{i}\})=\frac{1}{Z_{N}(\beta)}\,\exp[-\beta E[\{x_{i}\}]]$ can be identified with the joint distribution of $N$ real eigenvalues of an $N\times N$ matrix belonging to the Gaussian ensembles of the random matrix theory (RMT): respectively Gaussian Orthogonal Ensemble (GOE), Gaussian Unitary Ensemble (GUE) and Gaussian Symplectic Ensemble (GSE) Mehta_book ; Forrester_book . Gaussian ensembles are the cornerstones of RMT with myriads of applications, ranging from nuclear physics, mesoscopic transport, quantum chaos, number theory all the way to finance and big-data science Mehta_book ; Forrester_book ; Akemann_book ; LMV_book . The log-gas with arbitrary $\beta>0$ also appears in RMT as the joint PDF of the eigenvalues of the so called Dumitriu-Edelman $\beta$-ensemble of tridiagonal random matrices DE_2002 . The average density for large $N$ converges to the scaling form, independently of $\beta$,

[TABLE]

Thus the scaled density is supported over $[-\sqrt{2},+\sqrt{2}]$ and is known as the celebrated Wigner semi-circular law Wigner_1951 . This central result of RMT has been instrumental in understanding the global properties in a variety of systems including growth models in $1+1$ dimensions belonging to the KPZ universality class PS_2000 , non-interacting trapped fermions fermions_us1 ; fermions_us2 ; fermions_us3 where the Wigner semi-circle can also be obtained from the so-called local density approximation (LDA) LDA , non-intersecting Brownian motions fisher_vicious ; vicious ; BBMP_2014 , graph theory and communication networks Akemann_book .

The third physical example corresponds to $k=2$, i.e., with inverse-square repulsion. For $k=2$, (1) is the celebrated Calogero-Moser model Calogero (1969, 1971, 1975); Moser (1976) which is integrable and is ubiquitous in diverse fields Polychronakos (2006). Curiously, it turns out that for any finite $N$, in the minimum energy configuration of both the log-gas ($k\to 0^{+}$) and the Calogero-Moser model ($k=2$), the particle positions $x_{i}$’s coincide exactly Calogero_1981 ; AKD_2019 with the $N$ zeros of the Hermite polynomial of degree $N$. Consequently, for large $N$, the scaled average density for $k=2$ also approaches the Wigner semi-circular law in Eq. (3) and as in the log-gas case, the semi-circular law is independent of $\beta$. This is rather strange: even though the long-range repulsive interaction in the log-gas is much stronger than that of the inverse-square gas, the average density profile is identical in the two cases. This raises a very interesting and natural question: How does the shape of the average scaled density for large $N$ vary as one tunes the parameter $k$ ? The appearance of the semi-circular law both for $k\to 0^{+}$ and $k=2$ suggests an intriguing possibility of a ‘non-monotonic’ dependence of the shape of the average density as one tunes up $k$. This question on the dependence of the average density on $k$ also has practical implications in a number of physical contexts. For example, in a typical cold atom experimental set-up, a quadratic confining potential is natural due to the usage of optical laser traps. In addition, it is possible to induce long-range power-law repulsive interactions between such atoms. For instance, charged particles interacting via the $3d$ Coulomb repulsion, but confined to a line by highly anisotropic optical trap, would correspond to $k=1$. Similarly, $k=3$ describes a dipolar gas confined to $1d$ Lu et al. (2011); Griesmaier et al. (2005); Ni et al. (2010). For $k\to\infty$, the Riesz gas reduces to a truely short-ranged repulsive gas similar to the harmonically confined screened Coulomb or Yukawa-gas studied in Refs. Cunden_2017, ; Cunden_2018, . In addition, other values of $k$ have either been realized in experimental setups or could potentially be realized Brown and Carrington (2003).

In this Letter, we address this interesting question of the dependence of the average density on $k$ and obtain exact results for large $N$. We show that for large $N$ the average density is independent of $\beta$ for all $k$ (for sufficiently large $\beta$) footnote_beta , but has a rich and nontrivial dependence on $k$. On general grounds the average density is expected to have a scaling form $\langle\rho_{N}(x)\rangle\approx N^{-\alpha_{k}}\tilde{\rho}_{k}(x/N^{\alpha_{k}})$ for large $N$, where $N^{\alpha_{k}}$ corresponds to the typical scale of the position of the particles in the trap. Indeed, we do find this behavior, but with a twist. We show that there is a drastic change of behavior of the exponent $\alpha_{k}$ as well as the scaling function $\tilde{\rho}_{k}(y)$ at $k=1$. For the exponent, we get

[TABLE]

The scaling function $\tilde{\rho}_{k}(y)$ has a support over $[-\ell_{k}/2,+\ell_{k}/2]$ and can be expressed as $\tilde{\rho}_{k}(y)=\ell_{k}^{-1}F_{k}(y/\ell_{k})$ where $F_{k}(z)$ is given by

[TABLE]

with $-1/2\leq z\leq 1/2$ and $B(a,b)$ is the standard Beta function. Thus the density either diverges or vanishes at the two scaled edges $z=\pm 1/2$ with an exponent $\gamma_{k}$ which also exhibits a change of behavior at $k=1$, namely

[TABLE]

The support length $\ell_{k}$ is non-universal and depends explicitly on $k$ and the coupling strength $J$ (for the exact expressions of $\ell_{k}$, see Eqs. (36) and (65) of Supp. Mat. sup ). The scaling function $F_{k}(z)$ depends only on $k$, and is independent of $\beta$ and $J$. The case $k=1$ is marginal with additional logarithmic corrections (we discuss this later). We show that this change of behavior at $k=1$ can be traced back to the fact that, for $k<1$, the large distance behavior of the interaction term controls the large $N$ behavior of the density. In contrast, for $k>1$, the limiting density is determined by the short distance behavior of the interaction term. This gives rise to an effective field-theory that is fundamentally different for $k<1$ and $k>1$. Thus $k\to 0^{+}$ (log-gas) and $k=2$ (inverse-square gas) share the same average density profile, but the physics is rather different in the two cases. For $k\to-1$, we recover the flat density of the 1$d$OCP. Also, in the limit $k\to\infty$ we again get a flat density, consistent with the results for the $1d$ harmonically confined Yukawa gas Cunden_2017 ; Cunden_2018 . We also performed Monte-Carlo (MC) simulations for several values of $k$, finding excellent agreement with our analytical predictions (see Fig. 1).

Regime 1: $-2<k<1$: Assuming both terms in the energy (1) are of the same order for large $N$, the energy scale can be estimated as follows. Let the typical position of a particle scales as $x_{\rm typ}\sim N^{\alpha_{k}}$ for large $N$, where $\alpha_{k}$ is to be determined. Then the first term in (1) scales as $\sim N^{2\alpha_{k}+1}$, while the second term (where the double sum contains typically $N(N-1)\approx N^{2}$ terms) scales as $\sim J\,N^{2-k\,\alpha_{k}}$. Demanding they are of the same order fixes the exponent $\alpha_{k}=\frac{1}{k+2}$ (see the first line in Eq. (4)). Hence the total energy scales as $E\sim N^{2\alpha_{k}+1}\sim N^{\frac{4+k}{2+k}}$ in this regime. To find the configuration that dominates the partition function $Z_{N}(\beta)$ for large $N$, we generalise the method used for the log-gas ($k\to 0^{+}$ limit) Dyson_1962 ; DM_2006 ; DM_2008 ; Saff_book . It turns out to be convenient to express the coarse grained energy in terms of a macroscopic density $\rho_{N}(x)$ and use the relation $\sum_{i=1}^{N}f(x_{i})\approx N\,\int f(x)\,\rho_{N}(x)\,dx$, valid for any smooth function $f(x)$. Next, we rescale $x=N^{\alpha_{k}}\,y$ with $y\sim O(1)$. Under this rescaling, the density transforms as $\rho_{N}(x)\approx N^{-\alpha_{k}}{\tilde{\rho}_{k}}(y=x\,N^{-\alpha_{k}})$, where we assume $\tilde{\rho}_{k}(y)$ is smooth and normalizable, $\int{\tilde{\rho}_{k}}(y)dy=1$. Consequently, the coarse grained partition function for large $N$ can be expressed as a functional integral over the density field $\tilde{\rho}_{k}(y)$ (for details see Supp. Mat. sup )

[TABLE]

where the action $\Sigma[\tilde{\rho}_{k}(y)]$ is given by (see also Ref. Serfaty_book, for a rigorous proof in the cases $0<k<1$)

[TABLE]

Here $\mu$ is the Lagrange multiplier that enforces the constraint $\int\tilde{\rho}_{k}(y)dy=1$. Note that in the integrand of Eq. (7), we have only kept the leading order contributions to the energy. Both the entropy term (generated in going from microscopic configurations to the macroscopic density) as well as the short distance behavior of the interaction energy have been neglected, as they are of lower order in $N$ for $-2<k<1$. This is valid as long as $\beta\gg N^{-2\alpha_{k}}$ where $\alpha_{k}=1/(k+2)$ footnote_beta ; sup . Thus the effective action $\Sigma[\tilde{\rho}_{k}(y)]$ is manifestly non-local reflecting the long-range nature of the repulsive interaction. We will see later that this non-locality manifests only for $-2<k<1$. For large $N$, the partition function in Eq. (7) can then be evaluated by the saddle point method. Minimizing the action $\Sigma[\tilde{\rho}_{k}(y)]$ in (Harmonically confined particles with long-range repulsive interactions) with respect (w.r.t.) to $\tilde{\rho}_{k}(y)$ gives the saddle point equation for the optimal density

[TABLE]

This equation is valid over the support of $\tilde{\rho}_{k}(y)$. The density is clearly symmetric in $y$, hence the support is over $[-\ell_{k}/2,\ell_{k}/2]$ where $\ell_{k}$ is fixed using the normalization $\int_{-\ell_{k}/2}^{\ell_{k}/2}\rho_{k}(y)dy=1$. Taking a further derivative of (9) w.r.t. $y$ leads to a singular integral equation

[TABLE]

where $PV$ denotes the principal value which needs to be taken only for $k>0$. Also, for $k\to 0$, $J$ has to be rescaled such that $J\,|k|=1$. Solving this singular integral equation poses the main technical challenge. Fortunately, it turns out that for $-2<k<1$, this equation can be transformed into the well studied Sonin form  Sonin ; Popov_1982 ; Widom_1999 ; Buldyrev_Sonin ; Derrida_2007 ; Cividini_2017 ; Asaf_2019 , and can subsequently be inverted to obtain $\tilde{\rho}_{k}(y)$ explicitly sup . We then obtain the exact saddle point density $\tilde{\rho}_{k}(y)=\ell_{k}^{-1}F_{k}(y/\ell_{k})$ where $\ell_{k}$ is given in Supp. Mat. sup and $F_{k}(z)$ is given in Eq. (5) with $\gamma_{k}=(k+1)/2$. Note that the Sonin inversion formula also indicates that there is no physical solution (saddle point density) for $k>1$. Hence, this solution is valid only in the range $-2<k<1$. Furthermore, since $\beta$ appears only in the factor $\beta\,N^{(4+k)/(2+k)}$ outside the action $\Sigma[\tilde{\rho}_{k}(y)]$ in (7), it is clear that the saddle point density is independent of $\beta$: large $N$ is equivalent to large $\beta$. In addition, the average density $\langle\rho_{N}(x)\rangle$, for large $N$, clearly coincides with the saddle point density as the average over all possible densities is dominated by the saddle point. In Fig. 1 upper panels, we compare our theoretical predictions with MC simulations for three representative values of $k$ in the range $-2<k<1$ and find excellent agreement. Note that in the range $-2<k<-1$ the density diverges at the edges $\pm l_{k}/2$, while for $-1<k<1$ the density vanishes at the edges. Exactly at $k=-1$, the density is flat, consistent with the 1$d$OCP result.

Regime 2: $k\geq 1$: It turns out that, for $k>1$, the interaction term in (1) containing the double sum is dominated by particles which are very close to each other, i.e., almost nearest neighbours. As a result, the short distance properties of the interaction term plays a more dominant role compared with its long distance behavior. This leads to an effective field theory which is local in the density and is much simpler. To compute the effective coarse grained energy for large $N$, we then take a different path than the $k<1$ case (for details see Supp. Mat. sup ). First, it is convenient to order the particle positions so that $x_{i}$ increases with the label $i$ (this is fine since the energy (1) is invariant under permutation of labels). We then replace the discrete particle label $i$ by a continuous coordinate and the position $x_{i}$ is approximated by a smooth continuous function $x(i)$. Next, we approximate $x_{i}-x_{j}\approx(i-j)x^{\prime}(i)$ where $x^{\prime}(i)=dx(i)/di$ and we have kept only the first term in the Taylor expansion anticipating that it captures the leading short distance behavior. Our next step is to express $x^{\prime}(i)$ in terms of the local smooth macroscopic density $\rho_{N}(x)$ (normalized to unity). In fact, the local slope of the smooth function $x(i)$, i.e., $x^{\prime}(i)>0$ is simply the inverse of the number density $N\rho_{N}(x)$, i.e., $x^{\prime}(i)=1/(N\rho_{N}(x))$. Thus the double sum in (1) can be approximated, to leading order for large $N$, by $\sum_{i\neq j}|x_{i}-x_{j}|^{-k}=\sum_{i\neq j}|i-j|^{-k}[N\rho_{N}(x_{i})]^{k}$. The sum over $j$, for fixed $i$, is convergent for all $k>1$ and simply gives a factor $2\,\zeta(k)$, where $\zeta(k)=\sum_{n=1}^{\infty}n^{-k}$ is the Riemann zeta function. Furthermore, the sum over $i$ can be replaced by an integral using the relation $\sum_{i=1}^{N}f(x_{i})\approx N\int f(x)\rho_{N}(x)dx$ mentioned before. Using this relation in both terms of (1) leads to a coarse grained energy ${\cal E}\equiv{\cal E}[\rho_{N}(x)]$ sup

[TABLE]

which is completely local in the density $\rho_{N}(x)$, unlike (Harmonically confined particles with long-range repulsive interactions) for $-2<k<1$ that involved densities at two space-separated points. We then rescale $x\to x\,N^{-\alpha_{k}}$ and write $\rho_{N}(x)=N^{-\alpha_{k}}\tilde{\rho}_{k}(y=x\,N^{-\alpha_{k}})$. It is easy to see that for both terms in (11) to be of the same order in $N$ for large $N$, we need to choose $\alpha_{k}=k/(k+2)$, as stated in the second line of (4). Hence, the total energy scales as $E\sim N^{\frac{3k+2}{k+2}}$ for large $N$. The coarse grained partition function $Z_{N}(\beta)$ can then be written as

[TABLE]

where the action $\Sigma[\tilde{\rho}_{k}(y)]$ is given by

[TABLE]

with $\mu$ again denoting the Lagrange multiplier enforcing the normalization of the density. Note that we have kept the leading order contribution for large $N$ in the integrand in (12) and again neglected the entropy as well as subdominant singular terms which are of lower order in $N$. For large $N$, the integral (12) can again be evaluated by the saddle point method. Minimizing the action gives the saddle point equation

[TABLE]

Trivially solving this equation gives, $\tilde{\rho}_{k}(y)=(1/\ell_{k})\,F_{k}(y/\ell_{k})$, with support over $[-\ell_{k}/2,\ell_{k}/2]$ where $\ell_{k}=2\sqrt{2\mu}$ is fixed from the normalization and is given explicitly in Supp. Mat. sup . The scaling function $F_{k}(z)$ is then of the form in (5) with the exponent $\gamma_{k}=1/k$. The saddle point density coincides with the average density $\langle\rho_{N}(x)\rangle$ for large $N$. In addition, since $\beta$ appears only in the combination $\beta N^{(3k+2)/(k+2)}$ outside the action in (12), clearly the saddle point density and hence the average density $\langle\rho_{N}(x)\rangle$ is independent of $\beta$ for large $N$, as long as $\beta\gg N^{-2\alpha_{k}}$ where $\alpha_{k}=k/(k+2)$ footnote_beta ; sup . This analytical prediction is then verified in MC simulations (see the bottom panels in Fig. 1).

The marginal case $k=1$ lies at the borderline between Regime 1 and Regime $2$. In this case, one would expect logarithmic corrections. Indeed, we find sup that the average density approaches a scaling form, $\langle\rho_{N}(x)\rangle\sim L_{N}^{-1}\tilde{\rho}_{1}(x/L_{N})$, where the typical position of a particle scales as $L_{N}\approx(N\,\ln N)^{1/3}$ for large $N$. The scaling function $\tilde{\rho}_{1}(y)$ is supported over $[-\ell_{1}/2,\ell_{1}/2]$ with $\ell_{1}=2\sqrt{2\mu}$ and is given by

[TABLE]

This can be also cast in the scaling form $\tilde{\rho}_{1}(y)=(1/\ell_{1})\,F_{1}(y/\ell_{1})$ where $F_{1}(z)$ is given in (5) with $\gamma_{1}=1$. Numerical simulations are in good agreement with our analytical prediction, as shown in the Supp. Mat. sup .

Conclusions: In this Letter, we have computed analytically the average density profile of a classical gas of harmonically confined particles that repel each other with the repulsive interaction behaving as a power-law with exponent $-k$ of the distance between any pair of particles. Our result generalizes in a nontrivial way, to arbitrary $k>-2$, the three famous classical examples: the $1d$OCP ($k=-1$), the Dyson’s log-gas in RMT ($k\to 0^{+}$) and the Calogero-Moser model ($k=2$). We have shown that the underlying effective field theory that determines the average density profile for large $N$ is governed by fundamentally different physics for $-2<k<1$ and $k>1$. In the former case, the large distance behavior of the interaction potential dominates, while the latter case is governed by its short distance behavior. It would be interesting to study other observables beyond the average density for general $k>-2$. For instance, for the log-gas ($k\to 0^{+}$) the position of the rightmost particle $x_{\rm max}$ (the largest eigenvalue of a random matrix), centered and scaled, is known to converge to the celebrated Tracy-Widom distribution Tracy_Widom . The corresponding extreme value distribution for $k=-1$ has recently been computed exactly Dhar et al. (2017); dhar2018 and the case $k=2$ has been recently computed numerically AKD_2019 . It would be interesting to compute the limiting distribution of $x_{\rm max}$ for generic $k>-2$. Finally, it would be interesting to see if our predictions for the average density can be measured in cold atom experiments. From that perspective, it would be nice to extend our results for the density profile to higher dimensions.

Note added in proof: After submission of the work, we came to know from O. Zeitouni that the case $k>1$ was also studied recently in the mathematics literature HLSS18 .

Acknowledgments: We would like to thank T. Leblé, E. Saff and S. Serfaty for pointing out useful references and O. Zeitouni for stimulating discussions. MK would like to acknowledge support from the project 6004-1 of the Indo-French Centre for the Promotion of Advanced Research (IFCPAR), the Ramanujan Fellowship SB/S2/RJN-114/2016 and the SERB Early Career Research Award ECR/2018/002085 from the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India. AD, AK, SNM and GS would like to acknowledge support from the project 5604-2 of the Indo-French Centre for the Promotion of Advanced Research (IFCPAR). This work was supported by a research grant from the Center for Scientific Excellence at the Weizmann Institute of Science. AD, AK, SNM and GS acknowledge the hospitality of the Weizmann Institute during the SRITP Workshop “Correlations, Fluctuations and Anomalous Transport in Systems Far from Thermal Equilibrium” held at the Weizmann Institute in January 2018. SNM acknowledges the hospitality of the Weizmann Institute during a visit as a Weston Professor in 2019 and the support from the Science and Engineering Research Board (SERB, government of India), under the VAJRA faculty scheme (Ref. VJR/2017/000110) during a visit to Raman Research Institute, where part of this work was carried out. We would like to thank the ICTS program “Universality in random structures: Interfaces, Matrices, Sandpiles (Code: ICTS/URS2019/01)” for enabling valuable discussions with many participants. AK acknowledges support from DST grant under project No. ECR/2017/000634.