Finite temperature equilibrium density profiles of integrable systems in confining potentials

Jitendra Kethepalli [email protected] International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru – 560089, India Debarshee Bagchi [email protected] International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru – 560089, India Abhishek Dhar [email protected] International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru – 560089, India Manas Kulkarni [email protected] International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru – 560089, India Anupam Kundu [email protected] International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru – 560089, India

(December 4, 2024)

Abstract

We study the equilibrium density profile of particles in two one-dimensional classical integrable models, namely hard rods and the hyperbolic Calogero model, placed in confining potentials. For both of these models the inter-particle repulsion is strong enough to prevent particle trajectories from intersecting. We use field theoretic techniques to compute the density profile and their scaling with system size and temperature, and compare them with results from Monte-Carlo simulations. In both cases we find good agreement between the field theory and simulations. We also consider the case of the Toda model in which inter-particle repulsion is weak and particle trajectories can cross. In this case, we find that a field theoretic description is ill-suited due to the lack of a thermodynamic length scale. The density profiles for the Toda model obtained from Monte-Carlo simulations can be understood by studying the analytically tractable harmonic chain model (Hessian approximation of the Toda model). For the harmonic chain model one can derive an exact expression for the density that shines light on some of the qualitative features of the Toda model in a quadratic trap. Our work provides an analytical approach towards understanding the equilibrium properties for interacting integrable systems in confining traps.

I Introduction

Integrable classical systems [1, 2, 3] have a macroscopic number of constants of motion that are in involution with each other. In phase space, these systems (i) have regular periodic orbits (invariant torus), (ii) are characterized by zero Lyapunov exponents, and (iii) generally resist thermalization to a Gibbs state. However, many-body integrable systems are also believed to be extremely fragile in the presence of external perturbations, and become nonintegrable, ergodic, and chaotic, retaining only a few constants of motion [4]. Consequently, integrable systems are rare, and nonintegrability arising due to imperfections dominates the natural world. For example, most experiments [5, 6, 7] are performed in confining potentials where we expect that integrability will be lost and thermalization to occur. Recent theoretical studies have addressed thermalization and transport in such trapped integrable models [8, 9, 10, 11, 12, 13, 14, 15]. To study thermalization, one needs to have a clear understanding of the thermal equilibrium state. One simple characterization is to look at the equilibrium density profile of the particles in the trap which is the most commonly measured quantity in experiments [16, 17, 18].

In this work, we focus on equilibrium density profile of two one-dimensional short-range integrable classical models in the presence of integrability-breaking external potentials. The integrable models considered here are the gas of hard rods [19] and the hyperbolic Calogero model [20, 21, 22]. The external trap potential keeps the particles spatially confined and breaks integrability. Such systems have been studied recently and many surprising results have been reported. For example the gas of hard rods [11] and the Lieb-Liniger model [13] in quadratic trap were investigated. It was found that these systems do not thermalize even in the presence of the quadratic trap. Under out-of-equilibrium conditions, drastically slow relaxations to a nonequilibrium steady state and large finite size effects have also been observed for the Toda model with harmonic (quadratic) pinning potential [8, 12]. In another recent work [23], a similar observation was made for the nonlinearly perturbed Toda model and a universal scaling of the thermalization time has been reported. Studies of the integrable Calogero model in the presence of external confining potentials have also been undertaken in recent times, see for example Refs. [24, 10, 25, 26, 15].

The equilibrium properties of trapped interacting particles have recently been studied where field theoretic techniques are used to compute the equilibrium density profiles and fluctuations [27, 28, 29, 30, 31, 32, 33]. Here we adapt these field theoretic procedures to study the equilibrium properties of hard rods and the hyperbolic Calogero system in the presence of external trapping potentials. The field theory presented here predicts quite accurately the equilibrium density profile of these two models, and their scaling with system size and temperature, as obtained from Monte-Carlo (MC) simulations.

However, we find that the behavior of trapped integrable models with weak repulsion, such as the Toda model, appears to be strikingly different from the above-mentioned models. It turns out that for such short-range models where particle crossings are allowed, the density profile is localized on a length scale that is system size independent, thereby rendering the field theoretic description ill-suited. In a suitable parameter regime, the equilibrium density profiles of the Toda model can be understood using the Hessian approximation of the Toda interaction which is the nearest neighbor harmonic chain. Note that the harmonic chain model in the quadratic trap is analytically tractable and thereby provides a transparent way of understanding the density profiles.

The paper is organized as follows. We describe the models and definitions in Sec. II. Thereafter, in Sec. III, we present the field theory for the hard rods gas and the hyperbolic Calogero model in both quadratic and quartic traps. In Sec. IV we compute the densities and extract their scaling with system size and temperature for (i) hard rods gas in Sec. IV.1 and (ii) hyperbolic Calogero model in Sec. IV.2. We verify the analytical results using Monte-Carlo (MC) simulations. Next, in Sec. IV.3, we study the equilibrium density profiles of the Toda model along with its Hessian approximation in a suitable parameter regime. We summarize the main results in Sec. V and end with a discussion of open questions in such integrability broken classical systems. The Appendix is organised as follows. In Appendix A, we derive the field theory for (i) hard rods gas and (ii) hyperbolic Calogero model in external confining traps. In Appendix B, we compute the analytical form of the densities for low and high values of the temperature. In Appendix C, we derive the equilibrium density profile for the quadratically confined nearest neighbour harmonic chain.

II Models and definitions

We study two short-range models given by a Hamiltonian of the form

H(\{x_{i},p_{i}\})=\sum_{i=1}^{N}\left[\frac{p_{i}^{2}}{2m}+U_{\delta}(x_{i})\right]+\frac{1}{2}\sum_{i=1}^{N}\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}V(x_{i}-x_{j}),

(1)

where $\{x_{i},p_{i}\}$ are the position and momentum of the $i^{\rm th}$ particle ( $1\leq i\leq N$ ), each of mass $m$ which we set to unity. The second term on the right-hand side of Eq. (1) is the external potential

\displaystyle U_{\delta}(x)=\frac{x^{\delta}}{\delta},

(2)

which we take to be of quadratic ( $\delta=2$ ) or quartic ( $\delta=4$ ) form. The third term in Eq. (1) is the interaction term, which for hard rods (HR) of length $a$ is

V_{R}(r)=\begin{cases}0&\text{\leavevmode\nobreak\ for\leavevmode\nobreak\ }\quad r>a\\ \infty&\text{\leavevmode\nobreak\ for\leavevmode\nobreak\ }\quad r\leq a.\end{cases}

(3)

Note that in Eq. (3) the subscript ‘ $R$ ’ in $V_{R}(r)$ stands for the hard rods gas. For the hyperbolic Calogero (HC) model each particle is coupled to every other particles in the system with the interaction potential

V_{C}(r)=\frac{J}{\sinh^{2}|r|}.

(4)

In Eq. (4), the subscripts ‘ $C$ ’ in $V_{C}(r)$ stand for the hyperbolic Calogero model and $J>0$ is the strength of the repulsive interaction.

We consider these systems to be in their respective thermal equilibrium states described by the canonical Gibbs distribution

P(\{x_{i},p_{i}\})=\frac{e^{-\beta H(\{x_{i},p_{i}\})}}{Z_{\beta}(N)},

(5)

where $\beta=1/T$ is the inverse temperature and $Z_{\beta}(N)$ is the partition funtion. We are interested in the spatial density profile

\rho(x)=\sum_{i=1}^{N}\left<\delta(x-x_{i})\right>_{\beta},

(6)

where $\langle\ldots\rangle_{\beta}$ denotes the average over the thermal distribution given in Eq. (5). In particular, we will examine the dependence of the density profile $\rho(x)$ on system parameters, such as the number of particles $N$ and the temperature $T$ . In the following sections, we address these questions using field theory and MC simulations.

Refer to caption — Figure 1: A schematic representation for partitioning the system into $N_{b}$ subsystems, each of size $\Delta$ . Here $s=1,2,\ldots,N_{b}$ denotes the subsystem index. Note that the particles are represented by orange circles. In the $s^{\rm th}$ subsystem there are $n_{s}$ particles. The double arrow indicates the extent of each subsystem. While analyzing each subsystems we take large- $n_{s}$ and then for the complete system we finally take small- $\Delta$ limit.

III Field theory formalism

To obtain the thermal properties one needs to compute the partition function $Z_{\beta}(N)$ which is generally a hard task in microscopic variables. Therefore often one resorts to field theoretic (macroscopic) approach to compute $Z_{\beta}(N)$ . In this method, the partition function is written as a functional integral over density fields. This procedure has been commonly used in several contexts such as Landau theory [34], random matrix theory [27], general Coulomb gas [35], and long-range interacting particles [36, 29]. Despite this progress there has been only a few rigorous comparisons between densities and other equilibrium properties obtained from microscopic and macroscopic (field theory) computation [37, 29, 30, 31, 33, 32].

In this section, we describe a macroscopic procedure and construct a field theory adapted appropriately for our models. We start with the partition function

	$\displaystyle Z_{\beta}(N)$	$\displaystyle=\int_{-\infty}^{\infty}\prod_{i=1}^{N}dp_{i}\int_{-\infty}^{\infty}{\bf dx}_{N}{\rm exp}(-\beta H(\{x_{i},p_{i}\})),\leavevmode\nobreak\ \text{where}$
	$\displaystyle\int_{w}^{z}\leavevmode\nobreak\ {\bf dx}_{N}$	$\displaystyle\equiv\int_{w}^{z}dx_{1}\int_{x_{1}}^{z}dx_{2}\ldots\int_{x_{N-1}}^{z}dx_{N}.$		(7)

Since the position and momentum variables are uncoupled, the partition function reduces to

\displaystyle Z_{\beta}(N)=Z^{(\rm\mathcal{K})}_{\beta}(N)\leavevmode\nobreak\ Z^{(\rm\mathcal{C})}_{\beta}(N),

(8)

where the configurational contribution to the partition function is given by

	$\displaystyle Z^{(\rm\mathcal{C})}_{\beta}(N)=$	$\displaystyle\int_{-\infty}^{\infty}{\bf dx}_{N}\leavevmode\nobreak\ \times$
		$\displaystyle{\rm exp}\Big{(}-\beta\Big{[}\sum_{i=1}^{N}U_{\delta}(x_{i})+\frac{1}{2}\sum_{i=1}^{N}\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}V(x_{i}-x_{j})\Big{]}\Big{)},$		(9)

and contribution due to kinetic terms is

\displaystyle Z^{(\rm\mathcal{K})}_{\beta}(N)

\displaystyle=\left(\dfrac{2\pi}{\beta}\right)^{N/2}.

(10)

Performing the multiple integrals in Eq. (III) is a hard problem. However, for (i) short-range repulsive interactions that diverge at vanishing separation, and (ii) slowly varying confining potentials, one can approximate the full partition function as follows.

We divide the system into $N_{b}$ subsystems, as shown in Fig. 1, where each subsystem $s$ contains a large number of particles $n_{s}$ . Note that the size of each subsystem, denoted by $\Delta$ , is small enough compared to the actual size of the gas and large enough to contain many particles such that the change in potential energy between two successive boxes is smaller than thermal energy $T$ i.e., $|V_{\rm ext}(x_{s+1})-V_{\rm ext}(x_{s})|<T$ . The particles in each subsystem experience an effective constant potential that depends on the location of the subsystem $x_{s}$ inside the trap. The partition function $Z^{(\rm\mathcal{C})}_{\beta}(N)$ in Eq. (III) can be approximated (in the thermodynamic limit) as the product of the partition functions of these boxes,

\displaystyle Z^{(\rm\mathcal{C})}_{\beta}(N)

\displaystyle\approx{\rm exp}\left(\sum_{s=1}^{N_{b}}\log\big{[}Z_{\beta}(n_{s},x_{s},\Delta)\big{]}\right),

(11)

where the partition function of the $s^{\rm th}$ subsystem of size $\Delta$ , centered around $x_{s}$ containing $n_{s}$ particles is given by

	$\displaystyle Z_{\beta}(n_{s},x_{s},\Delta)=\int_{x_{s}-\frac{\Delta}{2}}^{x_{s}+\frac{\Delta}{2}}$	$\displaystyle{\bf dx}_{n_{s}}\prod_{i=1}^{n_{s}}\exp\Big{(}-\beta U_{\delta}(x_{i})\Big{)}$
		$\displaystyle\prod_{\begin{subarray}{c}i,j=1\\ j\neq i\end{subarray}}^{n_{s}}{\rm exp}\Bigg{(}-\beta\left[\frac{1}{2}V(x_{i}-x_{j})\right]\Bigg{)}$		(12)

The free energy per particle in the $s^{\rm th}$ box is given by

\displaystyle f\left(x_{s},\beta\right)=-\frac{1}{\beta n_{s}}\log\big{[}Z_{\beta}(n_{s},x_{s},\Delta)\big{]}.

(13)

We convert the summation in Eq. (11) over subsystem index $s$ to an integral over $x$ and get [see Appendix A]

\displaystyle\mathcal{F}[\rho(x),\beta]=\int_{-\infty}^{\infty}dx\leavevmode\nobreak\ \leavevmode\nobreak\ \rho(x)f\left(x,\beta\right).

(14)

where $\rho(x)$ is the density of particles at position $x$ . The free energy per particle $f(x_{s},\beta)$ defined in Eq. (13) can be computed from the partition function of the subsystem. As mentioned earlier, we assume that the subsystem size $\Delta$ is small enough such that all the $n_{s}$ particles with position $x_{i}$ (where $i=1,2..,n_{s}$ ), experience a constant potential $U_{\delta}(x_{i})\approx U_{\delta}(x_{s})$ . The subsystem partition function can then be approximated as

	$\displaystyle Z_{\beta}(n_{s},x_{s},\Delta)$	$\displaystyle\approx{\rm exp}\Big{(}-\beta n_{s}U_{\delta}(x_{s})\Big{)}\times$
		$\displaystyle\Bigg{[}\int_{x_{s}-\frac{\Delta}{2}}^{x_{s}+\frac{\Delta}{2}}{\bf dx}_{n_{s}}\prod_{\begin{subarray}{c}i,j=1\\ j\neq i\end{subarray}}^{n_{s}}{\rm exp}\Big{(}-\beta\left[\frac{1}{2}V(x_{i}-x_{j})\right]\Big{)}\Bigg{]}.$		(15)

Note that, in Eq. (III), the $x_{i}$ is a running integration variable not to be confused with the position of the center of the subsystem $x_{s}$ . The contribution to the free energy per particle from the $s^{\rm th}$ box is written as

\displaystyle f(x_{s},\beta)=U_{\delta}(x_{s})+f_{\rm int}(x_{s},\beta),

(16)

where

	$\displaystyle f_{\rm int}(x_{s},\beta)$	$\displaystyle=-\frac{1}{\beta n_{s}}\times$
		$\displaystyle\log\left(\int_{0}^{\Delta}{\bf dx}_{n_{s}}\prod_{\begin{subarray}{c}i,j=1\\ j\neq i\end{subarray}}^{n_{s}}{\rm exp}\Big{[}-\frac{\beta}{2}V(x_{i}-x_{j})\Big{]}\right).$		(17)

From Eq. (III) one can further rewrite $f_{\rm int}(x_{s},\beta)\equiv f_{\rm int}\big{(}\rho(x_{s}),\beta\big{)}$ . Furthermore, using Eq. (III) we can rewrite Eq. (14) as [see Appendix A]

\displaystyle\mathcal{F}[\rho,\beta]=\int_{-\infty}^{\infty}dx\leavevmode\nobreak\ \rho(x)\Big{\{}U_{\delta}(x)+f_{\rm int}\big{(}\rho(x),\beta\big{)}\Big{\}}.

(18)

For the HR and the HC models the explicit forms of the free energy are derived in Appendix A.1 and Appendix A.2 respectively.

The average thermal density can then be computed by extremizing the free energy in Eq. (18) with the constraint that the density is normalized

\displaystyle\int_{-\infty}^{\infty}dx\leavevmode\nobreak\ \rho(x)=N.

(19)

In the next section, we compute these densities for both HR and HC models and compare them with MC simulations.

IV Results from field theory and comparison with Monte-Carlo simulations

In this section, we adapt the field theory formalism discussed in Sec. III for the case of hard rods (HR) and hyperbolic Calogero (HC) model to compute the free energy. We extremize the obtained free energy along with the constraint that the density is normalized and this yields the average density profile.

IV.1 Hard rods model

For the HR model, the contribution due to interaction to the free energy per particle at position $x$ is given by (see Appendix A.1)

\displaystyle f_{\rm int}(\rho(x),\beta)=-\frac{1}{\beta}\log\left(\frac{1-a\leavevmode\nobreak\ \rho(x)}{\rho(x)}\right)+\frac{1}{\beta}.

(20)

Using Eq. (20) in Eq. (18) we get the free energy for the HR model [ignoring the density independent term $1/\beta$ in Eq. (20)]

\displaystyle\mathcal{F}_{R}\left[\rho(x),\beta\right]=\int_{-\infty}^{\infty}\leavevmode\nobreak\ dx\leavevmode\nobreak\ \rho(x)\Bigg{[}U_{\delta}(x)-\frac{1}{\beta}\log\left(\frac{1-a\leavevmode\nobreak\ \rho(x)}{\rho(x)}\right)\Bigg{]}.

(21)

The free energy in Eq. (21) is super-extensive i.e., $\mathcal{F}_{R}[\rho(x),\beta]\sim O(N^{\delta+1})$ since the ground state energy of $N$ hard rods in a confining potential $U_{\delta}(x)\sim x^{\delta}$ scales as $N^{\delta+1}$ . Therefore, the average thermal density $\rho^{*}(x,T)$ can be computed via saddle point approximation [29]. This amounts to extremizing the free energy along with the normalization constraint

\displaystyle\frac{\delta}{\delta\rho(x)}\mathcal{F}_{R}[\rho(x),\beta]\Bigg{|}_{\rho(x)=\rho^{*}(x,T)}=\mu_{N}(\beta),

(22)

where the chemical potential $\mu_{N}(\beta)$ is temperature dependent and can be extracted from normalization condition given in Eq. (19). Using Eq. (21) in Eq. (22) we get

	$\displaystyle\mu_{N}(\beta)=$	$\displaystyle U_{\delta}(x)$
		$\displaystyle-T\Bigg{[}\log\left(\frac{1-a\leavevmode\nobreak\ \rho^{}(x,T)}{\rho^{}(x,T)}\right)-\frac{1}{1-a\leavevmode\nobreak\ \rho^{*}(x,T)}\Bigg{]}.$		(23)

To obtain the system size dependence of the density profile, we define

\displaystyle\rho_{N}(x,T)=\frac{1}{N}\rho^{*}(x,T),

(24)

such that

\displaystyle\int_{-\infty}^{\infty}\leavevmode\nobreak\ dx\leavevmode\nobreak\ \rho_{N}(x,T)=1.

(25)

Using Eq. (24), Eq. (IV.1) can then be expressed as

	$\displaystyle\mu_{N}(\beta)=$	$\displaystyle U_{\delta}(x)$
		$\displaystyle-T\Bigg{[}\log\left(\frac{1-a\leavevmode\nobreak\ N\rho_{N}(x,T)}{N\rho_{N}(x,T)}\right)-\frac{1}{1-a\leavevmode\nobreak\ N\rho_{N}(x,T)}\Bigg{]}.$		(26)

To extract the system size ( $N$ ) and temperature ( $T$ ) dependence of the density $\rho_{N}(x,T)$ we substitute the following scaling form ansatz

\displaystyle\rho_{N}(x,T)

\displaystyle=N^{-\alpha_{R}}\leavevmode\nobreak\ \rho_{R}\left(y,c\right),\quad\leavevmode\nobreak\ \mu_{N}(\beta)=N^{\lambda_{R}}\leavevmode\nobreak\ \mu_{R}(c),

(27)

with the scaled variables given by

\displaystyle y

\displaystyle=\frac{x}{N^{\alpha_{R}}},\quad c=\frac{T}{N^{\gamma_{R}}},

(28)

in Eq. (IV.1). Here $\alpha_{R}$ and $\gamma_{R}$ are scaling exponents which are determined by requiring that Eq. (IV.1) is $N$ independent in the scaled variables. Doing so we get

\displaystyle\alpha_{R}=1,\quad\gamma_{R}=\delta\leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\quad\lambda_{R}=\delta.

(29)

The value $\alpha_{R}=1$ can be understood from the $O(N)$ extent of the density profile at zero temperature. This leads to $O(N^{\delta+1})$ energy of the system in the ground state. In order for the entropy term to contribute to the free energy one needs to scale the temperature by $N^{\delta}$ implying $\gamma_{R}=\delta$ . Eq. (IV.1) finally becomes

\displaystyle\mu_{R}(c)=\frac{y^{\delta}}{\delta}-c\Bigg{[}\log\left(\frac{1-a\leavevmode\nobreak\ \rho_{R}(y,c)}{\rho_{R}(y,c)}\right)-\frac{1}{1-a\leavevmode\nobreak\ \rho_{R}(y,c)}\Bigg{]}.

(30)

It is worth noting that the thermal equilibrium properties of hard rods in an external potential were studied in Ref. 38. Eq. (30) can be obtained from Eq. [13] of Ref. 38, when the density is assumed to vary slowly on the rod length scale $a$ . Since Eq. (30) is a transcendental equation, it is difficult to obtain an exact solution. We solve Eq. (30) numerically by fixing $\mu_{R}(c)$ such that the normalization constraint,

\displaystyle\int_{-\infty}^{\infty}dy\leavevmode\nobreak\ \rho_{R}(y,c)=1,

(31)

is satisfied. In Fig. 2 we show the comparison between the scaled density profile [obtained by solving Eq. (30)] and data from MC simulations (using the standard Metropolis algorithm) for three rescaled temperatures $c=0.1,1.0,10.0$ and three system sizes $N=32,64,128$ . We find quite remarkable scaling collapse of the MC data with system size which also agrees with the field theory results for both the quadratic ( $\delta=2$ ) and quartic ( $\delta=4$ ) traps.

Although explicit analytical solution of the saddle point given in Eq. (30) is highly nontrivial to obtain, one can study the behavior of the density for low $c\ll 1$ and high $c\gg 1$ analytically using asymptotic analysis (see Appendix B.1.1). At zero temperature the hard rods have a density profile given by

\displaystyle\rho_{R}(y,0)=\begin{cases}\frac{1}{a}&\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ y\leq\left|\frac{a}{2}\right|\\ 0&\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ y>\left|\frac{a}{2}\right|.\end{cases}

(32)

The density profile at low temperatures can then be approximated as

\displaystyle\rho_{R}(y,c)\stackrel{{\scriptstyle c\ll 1}}{{\approx}}\rho_{R}(y,0)+\rho_{1}(y,c),

(33)

where the deviation from the zero temperature density up to first iteration (see Appendix B.1.1 for more details) is given by

\displaystyle\rho_{1}(y,c)\approx\begin{cases}&-\frac{1}{a}\frac{c\leavevmode\nobreak\ \delta}{\left(y_{c}^{\delta}-y^{\delta}\right)}\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ |y|<y_{c}-O(c)\\ &\rho_{R}^{*}(1-a\leavevmode\nobreak\ \rho_{R}^{*})^{2}\times\\ &\left(1-\exp\left[-\frac{y_{c}^{\delta}-y^{\delta}}{c\leavevmode\nobreak\ \delta}\right]\right)\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ |y-y_{c}|<O(c)\\ &\frac{1}{e}\exp\left(\frac{y_{c}^{\delta}-y^{\delta}}{c\leavevmode\nobreak\ \delta}\right)\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ |y|>y_{c}+O(c).\end{cases}

(34)

Here $y_{c}$ is the position at which the term in the parenthesis of Eq. (30) changes sign and is given by

\displaystyle y_{c}=(\mu_{R}\delta)^{\frac{1}{\delta}}.

(35)

The density at $y=y_{c}$ is denoted by

\displaystyle\rho^{*}_{R}=\rho_{R}(y_{c},c).

(36)

Note that $y_{c}$ given in Eq. (35) determines the three regions (see Appendix B.1.1) (i) bulk region where $|y|<y_{c}-O(c)$ , (ii) edge region where $|y-y_{c}|\lesssim O(c)$ , and (iii) tail region where $|y|>y_{c}+O(c)$ , which are displayed in Eq. (34). The higher order corrections have also been computed and are presented in Appendix B.1.1. The expression Eq. (33) is verified with numerical solution of Eq. (30) for $c=0.01$ in Fig. 3 showing the three regions. For this comparison the value of the chemical potential $\mu_{R}(c)$ is taken from Fig. 4. Note that in Fig. 4 the behavior of $\mu_{R}(c)$ is non-monotonic: it increases initially as $c$ is increased from zero and thereafter decreases. This nonmonotonicity can be explained by noting that at smaller $c$ , particles can only be added to the edges of the system which requires more energy (owing to the confining potential), without a large increase in entropy. Hence $\mu_{R}(c)$ increases initially. However, for larger $c$ , the gas expands and this opens up gaps, larger than the size of the rods, between the particles in the bulk of the system. Consequently, one can easily add an extra rod with a small energy cost and a large entropy gain, essentially lowering the free energy change. Hence $\mu_{R}(c)$ decreases with $c$ for larger $c$ .

In the high temperature regime ( $c\gg 1$ ) the spread of the gas increases which in turn dilutes the gas i.e., $\rho_{R}(y,c)\ll 1$ . Using this low density approximation in Eq. (30), we obtain the approximate analytical expression of the density profile (up to the first iteration), given by (see Appendix B.1.2)

\displaystyle\rho_{R}(y,c)

\displaystyle\stackrel{{\scriptstyle c\gg 1}}{{\approx}}\frac{1}{e}\exp\left(\frac{\mu_{R}(c)}{c}-\frac{y^{\delta}}{c\delta}\right),

(37)

where the chemical potential $\mu_{R}(c)$ is obtained numerically by solving Eq. (30) along with the normalization condition [Eq. (31)] as shown in Fig. 4 and $e=2.71828$ is the Euler’s number. We can obtain higher order terms of the density by also considering subdominant corrections originating due to the presence of interaction as shown in Appendix B.1.2. The expression Eq. (37) and the subdominant corrections (up to third order) are verified with the numerical solution of Eq. (30) for both traps in Fig. 5 for $c=10.0$ .

IV.2 Hyperbolic Calogero model

Unlike the hard rods (HR) model, for the hyperbolic Calogero (HC) model, each particle is coupled to all the other particles. The field theoretic formulation of the hyperbolic Calogero model has been studied [26]. However, the average thermal density profiles at finite temperatures have not been computed yet. Based on the the approximate scheme outlined in Sec. III and the approach described in Refs. [39, 30, 29], we compute the finite temperature density profiles for the hyperbolic Calogero model below. The free energy in this case is given by (see Appendix A.2)

		$\displaystyle\mathcal{F}_{C}\left[\rho(x),\beta\right]$
		$\displaystyle=\int_{-\infty}^{\infty}\leavevmode\nobreak\ dx\leavevmode\nobreak\ \rho(x)\left[U_{\delta}(x)+J\zeta(2)\rho(x)^{2}+\frac{1}{\beta}\log\big{[}\rho(x)\big{]}\right],$		(38)

where $\zeta(k)=\sum_{n=1}^{\infty}n^{-k}$ is the Riemann Zeta function. Note that, despite the all-to-all coupling, the contribution to the free energy per particle due to interactions gets renormalized to a local term in the density field, and is given by

\displaystyle f_{\rm int}(\rho(x),\beta)=J\zeta(2)\rho(x)^{2}+\frac{1}{\beta}\log\rho(x).

(39)

Here $\beta^{-1}\log\rho(x)$ is the contribution due to the configurational entropy. To compute the average thermal density $\rho^{*}(x,T)$ , we extremize the free energy functional given in Eq. (IV.2) along with the normalization condition [Eq. (19)] which gives the chemical potential

\displaystyle\mu_{N}(\beta)=U_{\delta}(x)+3\zeta(2)\rho^{*}(x,T)^{2}+T\Big{(}1+\log\big{[}\rho^{*}(x,T)\big{]}\Big{)}.

(40)

As in the case of HR model, to obtain a scaling form for the density profile, we use the density normalized to unity $\rho_{N}(x,T)=\rho^{*}(x,T)/N$ in Eq. (40) and get

	$\displaystyle\mu_{N}(\beta)=U_{\delta}(x)$	$\displaystyle+3\zeta(2)N^{2}\rho_{N}(x,T)^{2}$
		$\displaystyle+T\Big{(}1+\log\big{[}N\rho_{N}(x,T)\big{]}\Big{)}.$		(41)

We can extract the system size ( $N$ ) and temperature ( $T$ ) dependence of the density $\rho_{N}(x,T)$ by substituting the scaling form ansatz

\displaystyle\rho_{N}(x,T)=N^{-\alpha_{C}}\rho_{C}\left(y,c\right),\quad\mu_{N}(\beta)=\mu_{C}(c)N^{\lambda_{C}},

(42)

with the scaled variables

\displaystyle y=\frac{x}{N^{\alpha_{C}}},\quad c=\frac{T}{N^{\gamma_{C}}},

(43)

in Eq. (IV.2). Here $\alpha_{C}$ and $\gamma_{C}$ are scaling exponents which are determined by requiring that Eq. (IV.2) is $N$ independent for large- $N$ and depends only on the scaling variables. Doing so, we get

\displaystyle\alpha_{C}=\frac{2}{2+\delta},\quad\gamma_{C}=\frac{2\delta}{2+\delta}\quad\text{and}\quad\lambda_{C}=\frac{2\delta}{2+\delta}.

(44)

The value $\alpha_{C}=2/(2+\delta)$ can be understood from the $O(N^{\frac{2}{2+\delta}})$ extent of the gas at zero temperature [29]. This leads to $O(N^{\alpha_{C}\delta+1})$ energy of the system in the ground state. In order for the entropy term to contribute to the free energy one needs to scale the temperature by $N^{\alpha_{C}\delta}$ implying $\gamma_{C}=\alpha_{C}\delta$ . Eq. (IV.2) finally becomes

\displaystyle\mu_{C}(c)=\frac{y^{\delta}}{\delta}+3\zeta(2)\rho_{C}(y,c)^{2}+c\log\rho_{C}(y,c).

(45)

We solve Eq. (45) numerically by fixing $\mu_{C}(c)$ such that the normalization constraint,

\displaystyle\int_{-\infty}^{\infty}dy\leavevmode\nobreak\ \rho_{C}(y,c)=1,

(46)

is satisfied. In Fig. 6, we show the comparison between the scaled density profile obtained by solving Eq. (45) and data from MC simulations for $c=0.1,1.0,10.0$ . We observe good agreement albeit with some small discrepancies, the origin of which is not understood clearly at present. The value of the chemical potential $\mu_{C}(c)$ is obtained as a function of $c$ by numerically solving Eq. (45) subject to the normalization condition Eq. (46), which is shown in Fig. 7. Unlike the HR model, we find that $\mu_{C}(c)$ decreases monotonically in this case.

Similar to the HR case, obtaining the exact solution of Eq. (45) is highly nontrivial for an arbitrary $c$ . However, we can study the average thermal density profiles using asymptotic analysis for small $c\ll 1$ and large $c\gg 1$ (see Appendix B.2). At zero temperature, $c=0$ , the density is exactly known and is governed by the interaction term only, since the contribution to the free energy from the entropy is zero. The zero temperature density is given by [29, 30]

\displaystyle\rho_{C}(y,0)=\begin{cases}A_{\delta}\left(l^{\delta}-y^{\delta}\right)^{\frac{1}{2}}&\text{for}\quad|y|<l\\ 0&\text{for}\quad|y|>l,\end{cases}

(47)

where

\displaystyle A_{\delta}=\left[3\delta\zeta(2)\right]^{-\frac{1}{2}}

(48)

and the edge of the support of the density is given by

\displaystyle l=\Big{(}\mu_{C}(0)\delta\Big{)}^{\frac{1}{\delta}}=\left(\frac{\delta}{2A_{\delta}{\rm B}\left(\frac{1}{\delta},\frac{3}{2}\right)}\right)^{\frac{2}{2+\delta}},

(49)

with

\displaystyle{\rm B}(x,y)=\int_{0}^{1}\leavevmode\nobreak\ dr\leavevmode\nobreak\ r^{x-1}(1-r)^{y-1},

(50)

being the Beta function. $\mu_{C}(0)$ in Eq. (49) is the scaled chemical potential at zero temperature, obtained by imposing the normalization condition [Eq. (46)], and is given by

\displaystyle\mu_{C}(0)=\frac{l^{\delta}}{\delta}=\left(\frac{\pi}{2}\right)^{\frac{\delta}{\delta+2}}\left(\frac{\delta^{-1/\delta}\Gamma\left(\frac{3}{2}+\frac{1}{\delta}\right)}{\Gamma\left(1+\frac{1}{\delta}\right)}\right)^{\frac{2\delta}{\delta+2}},

(51)

where

\displaystyle\Gamma[n]=\int_{0}^{\infty}\leavevmode\nobreak\ dx\leavevmode\nobreak\ x^{n-1}e^{-x},

(52)

is the Gamma function.

For $c\neq 0$ , the entropy starts contributing to the density. As the rescaled temperature is increased from zero, i.e., $c\ll 1$ , we can obtain the approximate analytical form of the density profile, as shown in the Appendix B.2.1, which is given by

\displaystyle\rho_{C}(y,c)\stackrel{{\scriptstyle c\ll 1}}{{\approx}}\rho_{C}(y,0)+\rho_{1}(y,c).

(53)

Here the deviation from zero temperature density (up to first iteration) is given by

\displaystyle\rho_{1}(y,c)\approx\begin{cases}\rho_{C}(y,0)\times&\\ \frac{\mu_{C}(c)-\mu_{C}(0)-c\log\rho_{C}(y,0)}{c+6\zeta(2)\rho_{C}(y,0)^{2}},&\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ |y|<y_{c}-O(c),\\ \frac{\rho_{C}^{*}\left(y_{c}^{\delta}-y^{\delta}\right)}{\delta\left(c+6\zeta(2)\rho_{C}^{*2}\right)}&\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ |y-y_{c}|<O(c),\\ \exp\left(\frac{y_{c}^{\delta}-y^{\delta}}{c\delta}\right)&\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ |y|>y_{c}+O(c),\end{cases}

(54)

and the higher order corrections are provided in Appendix B.2.1. Similar to the HR model, here $y_{c}=(\mu_{C}(c)\delta)^{\frac{1}{\delta}}$ and $\rho_{C}^{*}=\rho_{C}(y_{c},c)$ is the value of the density at $y=y_{c}$ . In Fig. 8, we find a good agreement between the expression Eq. (53) and the numerical solution of Eq. (45) for $c=0.01$ . Note that, for this comparison the value of the chemical potential $\mu_{C}(c)$ is taken from Fig. 7, where we recall that $\mu_{C}(c)$ is obtained by solving Eq. (45) along with the normalization condition [Eq. (46)].

As temperature increases the particles spread spatially over a wider region. Therefore, at high temperatures, $c\gg 1$ , the gas becomes dilute i.e. $\rho_{C}(y,c)\ll 1$ . Using this low density approximation in Eq. (45) yields (see Appendix B.2.2)

\displaystyle\rho_{C}(y,c)\stackrel{{\scriptstyle c\gg 1}}{{\approx}}\exp\left(\frac{\mu_{C}(c)}{c}-\frac{y^{\delta}}{c\delta}\right),

(55)

where $\mu_{C}(c)$ is obtained numerically from Fig. 7. The form of the density in Eq. (55) comes from the entropy which provides the dominant contribution to the density for $c\gg 1$ . In Fig. 9, for $c=10$ , we find a good agreement between the approximate expression of the density profile given in Eq. (55) (see Appendix B.2.2 for higher correction) and the numerical solution of the Eq. (45).

IV.3 Integrable models with particle-crossing

In the previous sections, we studied hard rods (HR) and hyperbolic Calogero (HC) models which have strong interparticle repulsions that prevent their trajectories from crossing. In this section, we study two models, namely Toda and Harmonic chains, which allow for crossing of particle trajectories due to their weak interparticle repulsion. For these models, the interactions are nearest neighbor. The Toda model takes the form

\displaystyle V_{T}(r_{i})=J\exp\left(-\frac{r_{i}}{d}\right),

(56)

where $r_{i}=x_{i+1}-x_{i}$ , $J>0$ (we set $J=1$ ) is the interaction strength and $d$ determines the length scale of the interaction. The subscript $T$ in the Eq. (56) stands for the Toda model. For the case of harmonic chain, the interaction takes the form

\displaystyle V_{H}(r_{i})=\frac{k_{2}}{2}r_{i}^{2}-k_{3}r_{i},

(57)

where $k_{2}$ and $k_{3}$ are the interaction strengths. The subscript $H$ in Eq. (57) stands for the harmonic chain. Note that in the absence of the trap these models are integrable, similar to HR and HC. Furthermore note that the form of the harmonic chain interaction in Eq. (57) is obtained by expanding the Toda interaction given in Eq. (56) up to the quadratic order in the nearest neighbor separation $r_{i}$ , with

\displaystyle k_{2}=\frac{1}{d^{2}}\quad\text{and}\quad k_{3}=\frac{1}{d}.

(58)

This is an analytically tractable model even in the presence of a quadratic trap given by $U_{\delta}(x)=k_{1}x^{\delta}/\delta$ with $\delta=2$ . We set $k_{1}=1$ .

In contrast to the HR [Sec. IV.1] and HC [Se . IV.2] models, we find that the field theoretic description fails to describe the equilibrium properties of the trapped Toda and harmonic chain models. The failure of the field theory in this case can be ascribed to the fact that the particles stay confined to a region of length $\sim O(N^{0})$ or smaller, due to the lack of strong repulsion in presence of external confining trap. To understand this, we consider the behavior of the system at zero temperature. By minimizing the energy one can find the particle positions and it turns out that for both the models a large number [ $\sim O(N)$ ] of particles are confined to a distance of $\sim O(N^{0})$ around the minimum of the external potential. This fact can be understood as follows. Assuming there is a length scale of order $O(N^{\alpha})$ , i.e., the particle positions are of order $x_{i}\sim N^{\alpha}y_{i}$ , we compute the contributions from the potential and interaction energy. The contribution from the external potential is of order $\sum_{i=1}^{N}U_{\delta}(x_{i})\sim O(N^{\alpha\delta+1})$ . The contribution from the interaction energy for the Toda model scales as $\sum_{i=1}^{N-1}V_{T}(x_{i+1}-x_{i})\sim O\big{(}N\exp\left({-N^{\alpha}}\right)\big{)}$ and for the Harmonic chain it scales as $\sum_{i=1}^{N-1}V_{H}(x_{i+1}-x_{i})\sim O(N^{2\alpha+1})$ . Comparing these energy contributions, we obtain $\alpha=0$ , for both models and for all values of $\delta>0$ , implying a length scale of $O(N^{0})$ . As a consequence the total energy is of $O(N)$ .

As the temperature $T$ is increased, maintaining $T\sim O(1)$ , the particles still stay extended over a region $\sim O(N^{0})$ . This is because the contribution to the free energy due to entropy is $O(N)$ which is of the same order as the energy contribution. This is in sharp contrast to the HR [Sec. IV.1] and HC [Sec. IV.2] models as discussed in the previous sections where the length scale increases with $N$ . Hence a field theory construction does not make sense for both the Toda model and the Harmonic chain due to the lack of a macroscopic length scale.

Therefore in order to understand the equilibrium properties of such models, we have performed detailed MC simulations of the Toda model with quadratic and quartic traps. In Fig. 10 we plot the density profile for the Toda model at $T=1$ for $d=0.01,1.0,10.0$ with $N=32,64,128$ . For $d=0.01,1.0,10.0$ , we find that the density profiles converge with increasing system size and have $O(N^{0})$ spread as expected. However, for $d=0.1$ we do not observe the convergence for the system sizes studied.

It is interesting to consider two special limits in the interaction length scale ( $d$ ) in the Toda model. For small $d\ll 1$ , the Toda interactions become similar to hard point gas, and for large $d\gg 1$ it can be approximated by nearest neighbor weakly interacting harmonic chain. In these two limits, the density profile has the form $\sim e^{-\beta U_{\delta}(x)}$ . More precisely in our case we get

\displaystyle\rho_{T}(x,\beta)\approx\frac{\beta^{\frac{1}{\delta}}}{\delta^{\frac{1}{\delta}}\Gamma\left[1+1/\delta\right]}{\rm exp}\left(-\beta\frac{x^{\delta}}{\delta}\right),

(59)

for both quadratic $(\delta=2)$ and quartic trap $(\delta=4)$ as verified in Fig. 11. Furthermore, note that the free energy of the Toda model $\mathcal{F}_{T}(\beta,d,J)$ at temperature $1/\beta$ can be related to the free energy of a Toda model at $\beta=1$ as

\displaystyle\mathcal{F}_{T}(\beta,d,J)=\mathcal{F}_{T}(1,d\beta^{\frac{1}{\delta}},\beta J)+\frac{N}{\beta\delta}\ln\beta.

(60)

Interestingly, Eq. (60) implies that the Toda model at any temperature can be mapped to Toda model at temperature $1/\beta=1$ with rescaled interaction strength $J\to\beta J$ and interaction length scale $d\to d\beta^{\frac{1}{\delta}}$ . Therefore studying density profiles for various values of $d$ is equivalent to studying the density profiles for different values of temperatures. Note that, the Toda model at high temperatures ( $\beta\ll 1$ ) can be approximated as a hard point gas ( $d\ll 1$ ) and for low temperatures ( $\beta\gg 1$ ) it can be approximated as a harmonic chain with nearest neighbor interactions ( $d\gg 1$ ) under an external pressure.

To better understand the features of the density profiles of the Toda model we study its harmonic limit [Eq. (57) with $k_{2}=1/d^{2}$ and $k_{3}=1/d$ ], which is analytically tractable for the quadratic trap. We obtain the exact expression (see Appendix C) for the density profile of the harmonic chain in a quadratic trap which is given by (see Appendix C)

\displaystyle\rho_{H}\left(x,\beta\right)

\displaystyle=\frac{1}{N}\sum_{l=1}^{N}\frac{1}{\sqrt{2\pi{\rm Var}(x_{l})}}{\rm exp}\left(-\frac{\Big{(}x-\langle x_{l}\rangle_{\beta}\Big{)}^{2}}{2\leavevmode\nobreak\ {\rm Var}(x_{l})}\right),

(61)

where the mean position of the $l^{\rm th}$ particle $\langle x_{l}\rangle_{\beta}$ and its variance ${\rm Var}(x_{l})$ are given in Eq. (169) and Eq. (170) of Appendix C respectively. This is plotted in Fig. 12 and compared with the corresponding MC simulation of the Toda model for $d=1.0$ and $10.0$ . We find good agreements both for $d=1.0,10.0$ . Note that already for $N=32,64,128$ the density profiles have converged to an $N$ -independent form. As mentioned earlier for the Toda model at $d=0.1$ , a slower convergence with $N$ was observed for the density profile (see Fig 10 b, f). Interestingly a similar slow convergence can be analytically demonstrated (see Appendix C) for a stiff (large $k_{2}$ ) harmonic chain. In Fig. 13, the density profiles for the harmonic chain with $T=1$ , $k_{1}=1$ , $k_{2}=100$ and $k_{3}=10$ for $N=32,64,128,256$ and $N\to\infty$ are shown. It can be seen that for increasing $N$ the density profiles converges slowly to the $N\to\infty$ curve.

Model \ Trap	Harmonic ( $\delta=2$ )	Quartic ( $\delta=4$ )
Hard rods	Eq. (27),	Eq. (27),
	$\alpha_{R}=1$ , $\gamma_{R}=2$ ,	$\alpha_{R}=1$ , $\gamma_{R}=4$ ,
	$\lambda_{R}=2$	$\lambda_{R}=4$
	Figs. 2a, 2b, 2c	Figs. 2d, 2e, 2f
Hyperbolic Calogero	Eq. (42),	Eq. (42),
	$\alpha_{C}=\frac{1}{2}$ , $\gamma_{C}=1$ ,	$\alpha_{C}=\frac{1}{3}$ , $\gamma_{C}=\frac{4}{3}$ ,
	$\lambda_{C}=1$	$\lambda_{C}=\frac{4}{3}$
	Figs. 6a, 6b, 6c	Figs. 6d, 6e, 6f

Table 1: A summary of the scaling behavior of the densities for the hard rods (HR) and the hyperbolic Calogero (HC) model in quadratic and quartic traps.

V Conclusion

To summarize, we have presented the equilibrium density profiles at finite temperatures of two integrable models, the hard rods and the hyperbolic Calogero model, in quadratic and quartic traps. For these models inter-particle repulsion is strong enough to avoid particle trajectories from crossing. The trap confines these systems spatially and breaks integrability. For these two models, we studied equilibrium density profiles using a field theory approach and Monte-Carlo (MC) simulations.

We developed appropriate field theory for these two models by extending the approaches used in Ref. [29]. From the field theory we computed the equilibrium density profiles, and their dependence on system size $N$ and temperature $T$ . The field theory calculations predict precise scaling forms for the equilibrium density profiles with respect to $N$ and $T$ . A summary of the scaling forms are given in Table 1. We find that the predictions from field theory for hard rods agree remarkably well with MC simulations (Fig. 2). For the hyperbolic Calogero model the agreement is also reasonably good (Fig. 6). On the other hand for integrable models that allow crossing of particle trajectories, such as the Toda model in quadratic and quartic trap, a field theoretic description is ill-suited due to the lack of a macroscopic length scale. For this case, we have presented microscopic analytical calculations, by employing Hessian approximation, and results from MC simulations.

Our work provides a framework for investigating the non-equilibrium dynamics, thermalization and transport in integrable models confined in external potentials. More precisely, one can ask whether these systems under Hamiltonian dynamics are ergodic, chaotic, and whether or not they equilibrate/thermalize, when placed in different confining traps. This is an area of active current research both theoretically [15] and experimentally [5].

Acknowledgements

M.K. would like to acknowledge support from the project 6004-1 of the Indo-French Centre for the Promotion of Advanced Research (IFCPAR), Ramanujan Fellowship (SB/S2/RJN-114/2016), SERB Early Career Research Award (ECR/2018/002085) and SERB Matrics Grant (MTR/2019/001101) from the Science and Engineering Research Board (SERB), Department of Science and Technology (DST), Government of India. A.K. acknowledges the support of the core research grant CRG/2021/002455 and the MATRICS grant MTR/2021/000350 from the SERB, DST, Government of India. A.D., M.K., and A.K. acknowledges support of the Department of Atomic Energy, Government of India, under Project No. 19P1112R&D.

Appendix A Derivation of the free energy

To obtain the free energy $\mathcal{F}[\rho,\beta]$ in Eq. (18), one first needs to compute $f_{\rm int}(\rho(x),\beta)$ defined in Eq. (III), where we recall that $f_{\rm int}(\rho(x),\beta)$ is the contribution to the free energy of the subsystem (recall Fig. 1) due to interactions (i.e., excluding the external confining potential). In the following, we present the calculation of $f_{\rm int}(\rho(x),\beta)$ for the hard rods (HR) model in Appendix A.1 and the hyperbolic Calogero (HC) model in Appendix A.2 separately.

A.1 Free energy for hard rods

The free energy per particle for hard rods of length $a$ , $f_{\rm int}(\rho(x),\beta)$ , can be calculated using the partition function [term in the parenthesis (square bracket) of Eq. (III)]

		$\displaystyle Z_{\rm int}(n_{s},x_{s},\Delta,\beta)$
		$\displaystyle=\int_{x_{s}-\Delta/2+\frac{a}{2}}^{x_{s}+\Delta/2-(n_{s}-\frac{1}{2})a}dy_{1}...\int^{x_{s}+\Delta/2-(n_{s}-i+\frac{1}{2})a}_{y_{i-1}+a}dy_{i}...\times$
		$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\int^{x_{s}+\Delta/2-\frac{a}{2}}_{y_{n_{s}-1}+a}dy_{n_{s}}$
		$\displaystyle=\int_{\frac{a}{2}}^{\Delta-(n_{s}-\frac{1}{2})a}dy_{1}...\int^{\Delta-(n_{s}-i+\frac{1}{2})a}_{y_{i-1}+a}dy_{i}...\int^{\Delta-\frac{a}{2}}_{y_{n_{s}-1}+a}dy_{n_{s}},$		(62)

where $y_{i}$ is the position of the $i^{\rm th}$ rod of the subsystem which is centered at $x_{s}$ and has a size $\Delta$ . In each subsystem there are $n_{s}$ hard rods. Note that, since the integrand in the second line of Eq. (A.1) is constant and translationally invariant, we have shifted the limits of the integrals from $y_{i}\to y_{i}-(x_{s}-\Delta/2)$ . This shift results in the integrals given in the third line of Eq. (A.1). These integrals can be computed using the variable transformation $z_{i}=y_{i}-(i-\frac{1}{2})a$ , which gives

\displaystyle Z_{\rm int}(n_{s},x_{s},\Delta,\beta)

\displaystyle={\rm exp}\Bigg{[}n_{s}\log\left(\frac{1-\rho(x_{s})a}{\rho(x_{s})}\right)-n_{s}\Bigg{]},

(63)

where we introduce the density in the given subsystem

\displaystyle\rho(x_{s})=\frac{n_{s}}{\Delta}.

(64)

The free energy per particle in a given subsystem in the large $n_{s}$ limit is given by [19]

	$\displaystyle f_{\rm int}\left(x_{s},\beta\right)$	$\displaystyle=-\frac{1}{n_{s}\beta}\log\Big{[}Z_{\rm int}(n_{s},x_{s},\Delta,\beta)\Big{]}$
		$\displaystyle=-\frac{1}{\beta}\log\left[\frac{1-a\rho(x_{s})}{\rho(x_{s})}\right]+\frac{1}{\beta}.$		(65)

One can see that, from the partition function in Eq. (A.1), the logarithmic term in Eq. (A.1) is the configurational entropy which includes the effect of hard rod exclusion. Note that the free energy due to interaction is a function of the density field and we rewrite the arguments of

\displaystyle f_{\rm int}(x_{s},\beta)

\displaystyle\equiv f_{\rm int}(\rho(x_{s}),\beta).

(66)

The total (i.e. including the contribution due to the external potential) free energy of the entire system $\mathcal{F}_{R}[\rho(x_{s}),\beta]$ is obtained by summing over the total free energy

\displaystyle n_{s}f(x_{s},\beta)=n_{s}f_{\rm int}(\rho(x_{s}),\beta)+n_{s}U_{\delta}(x_{s}),

(67)

associated with each subsystem. We therefore get

	$\displaystyle\mathcal{F}_{R}[\rho(x_{s}),\beta]$	$\displaystyle=\sum_{s=1}^{N_{b}}n_{s}f(x_{s},\beta),$
		$\displaystyle=\sum_{s=1}^{N_{b}}\leavevmode\nobreak\ \rho(x_{s})\leavevmode\nobreak\ \Delta\leavevmode\nobreak\ \Big{[}f_{\rm int}(\rho(x_{s}),\beta)+U_{\delta}(x_{s})\Big{]}.$		(68)

Converting the summation to integration i.e., $\sum_{s=1}^{N_{b}}\Delta\to\int_{-\infty}^{\infty}dx$ we obtain

\displaystyle\mathcal{F}_{R}[\rho(x),\beta]

\displaystyle=\int_{-\infty}^{\infty}dx\leavevmode\nobreak\ \rho(x)\Big{[}f_{\rm int}(\rho(x),\beta)+U_{\delta}(x)\Big{]}

(69)

Using Eq. (A.1) in Eq. (69), we obtain

\displaystyle\mathcal{F}_{R}\left[\rho(x),\beta\right]=\int_{-\infty}^{\infty}dx\leavevmode\nobreak\ \rho(x)\Bigg{[}U_{\delta}(x)-\frac{1}{\beta}\log\left(\frac{1-a\rho(x)}{\rho(x)}\right)\Bigg{]},

(70)

which is the free energy of the hard rods in an external trap $U_{\delta}(x)$ given in Eq. (21) of the main text. In Eq. (70) we have ignored the density independent term from Eq. (A.1).

A.2 Free energy for hyperbolic Calogero model

The field theoretic description of the hyperbolic Calogero model in external traps is well understood [26]. In this section we present an alternative derivation of the total free energy $\mathcal{F}_{C}[\rho(x),\beta]$ of the system. Using the approximate scheme described in Refs. [39, 29] we compute the free energy per particle $f_{\rm int}(\rho(x),\beta)$ of the subsystem due to the interaction which is described below.

The free energy per particle for the HC model, $f_{\rm int}(\rho(x),\beta)$ , can be calculated using the partition function which we recall to be

	$\displaystyle Z_{\beta}(n_{s},x_{s},\Delta)$	$\displaystyle\approx{\rm exp}\big{(}-\beta n_{s}U_{\delta}(x_{s})\big{)}\times$
		$\displaystyle\Bigg{[}\int_{x_{s}-\frac{\Delta}{2}}^{x_{s}+\frac{\Delta}{2}}{\bf dx}_{n_{s}}\prod_{\begin{subarray}{c}i,j=1\\ j\neq i\end{subarray}}^{n_{s}}{\rm exp}\Big{(}-\frac{\beta}{2}\left[V(x_{i}-x_{j})\right]\Big{)}\Bigg{]}.$		(71)

For the HC model Eq. (A.2) becomes

\displaystyle Z_{\beta}(n_{s},x_{s},\Delta)

\displaystyle\approx{\rm exp}\big{(}-\beta n_{s}U_{\delta}(x_{s})\big{)}Z_{\rm int}(n_{s},x_{s},\Delta,\beta),

where

	$\displaystyle Z_{\rm int}(n_{s},x_{s},\Delta,\beta)=\int_{0}^{\Delta}$	$\displaystyle{\bf dx}_{n_{s}}\times$
		$\displaystyle{\rm exp}\left[-\frac{\beta J}{2}\sum_{i=1}^{n_{s}}\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{n_{s}}\frac{1}{\sinh^{2}\left(\|x_{i}-x_{j}\|\right)}\right],$		(72)

where $x_{i}$ is the position of the $i^{\rm th}$ particle and $J$ is the interaction strength. As mentioned in the main text [Sec. III], the $x_{i}$ is a running integration variable not to be confused with the position of the center of the subsystem $x_{s}$ . One can approximate the exponential term in the integrand of Eq. (A.2) as

	$\displaystyle{\rm exp}\left(-\frac{\beta J}{2}\sum_{i=1}^{n_{s}}\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{n_{s}}\frac{1}{\sinh^{2}\left(\|x_{i}-x_{j}\|\right)}\right)$
	$\displaystyle\approx{\rm exp}\left(-\frac{\beta J}{2}\frac{n_{s}^{2}}{\Delta^{2}}\sum_{i=1}^{n_{s}}\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{n_{s}}\frac{1}{\left(\|i-j\|\right)^{2}}\right),$
	$\displaystyle\approx{\rm exp}\left(-\beta J\frac{n_{s}^{3}}{\Delta^{3}}\zeta(2)\Delta\right),$		(73)

where in the second line of Eq. (73) we approximated $x_{i}\approx i\Delta/n_{s}$ for all $i$ since $\Delta$ is assumed to be small enough to ensure uniform density over the subsystem. We have also used $\zeta(2)=\sum_{i=1}^{\infty}1/i^{2}$ , where $\zeta(k)=\sum_{i=1}^{\infty}1/i^{k}$ represents the Riemann zeta function. Using Eq. (73) the partition function in Eq. (A.2) takes the form

	$\displaystyle Z_{\rm int}(n_{s},x_{s},\Delta,\beta)$	$\displaystyle\approx{\rm exp}\left(-\beta J\frac{n_{s}^{3}}{\Delta^{3}}\zeta(2)\Delta\right)\int_{0}^{\Delta}{\bf dx}_{n_{s}},$
		$\displaystyle={\rm exp}\left(-\beta J\frac{n_{s}^{3}}{\Delta^{3}}\zeta(2)\Delta\right)\frac{\Delta^{n_{s}}}{n_{s}!},$		(74)

which can be rewritten using Stirling’s approximation $\log[n!]\approx n\leavevmode\nobreak\ \log[n]-n$ as

	$\displaystyle Z_{\rm int}$	$\displaystyle(n_{s},x_{s},\Delta,\beta)$
		$\displaystyle\asymp{\rm exp}\Bigg{(}-n_{s}\Big{[}\log\big{[}\rho(x_{s})\big{]}+\zeta(2)\beta J\rho(x_{s})^{2}\Big{]}\Bigg{)},$		(75)

where we recall that $\rho(x_{s})=n_{s}/\Delta$ . Hence, using the first line in Eq. (A.1), the free energy per particle of the subsystem becomes

\displaystyle f_{\rm int}\left(\rho(x_{s}),\beta\right)=J\zeta(2)\rho(x_{s})^{2}+\frac{1}{\beta}\log\left[\rho(x_{s})\right].

(76)

Similar to procedure detailed in Appendix A.1, using the above expression Eq. (76) we can compute the total free energy of the system as

	$\displaystyle\mathcal{F}_{C}\left[\rho(x),\beta\right]=\int_{-\infty}^{\infty}\leavevmode\nobreak\ dx\leavevmode\nobreak\ \rho(x)$	$\displaystyle\Bigg{(}U_{\delta}(x)+J\zeta(2)\rho(x)^{2}$
		$\displaystyle+\frac{1}{\beta}\log\rho(x)\Bigg{)},$		(77)

which is the expression for the free energy [see Eq. (IV.2) of the main text] of the HC model in an external trap $U_{\delta}(x)$ .

Appendix B Analytical forms of density profiles for hard rods and hyperbolic Calogero model at low and high rescaled temperatures $c$

To obtain the exact analytical expression for the equilibrium density profiles by solving the transcendental equations, Eq. (30) for HR model and Eq. (45) for HC model, is a highly non-trivial. However, we can obtain approximate expressions for the densities at low ( $c\ll 1$ ) and high ( $c\gg 1$ ) rescaled temperatures $c$ . For $c\ll 1$ , this is done by approximating the equilibrium density profile to be a small deviation around the zero temperature density profile. On the other hand, for $c\gg 1$ , the particles expected to spread far apart, thereby diluting the gas. Thus in this regime ( $c\gg 1$ ) it is reasonable to assume the density to be very small. In this section, using the above mentioned assumptions for the saddle point equations i.e., Eq. (30) for HR model and Eq. (45) for HC model, we discuss the analytical forms of the density profiles.

B.1 Analytical forms of the density profiles for hard rods

In this subsection, we discuss the case of hard rods, recall that the saddle point equation for the hard rods is

\displaystyle\mu_{R}(c)=\frac{y^{\delta}}{\delta}-c\Bigg{[}\log\left(\frac{1-a\leavevmode\nobreak\ \rho_{R}(y,c)}{\rho_{R}(y,c)}\right)-\frac{1}{1-a\leavevmode\nobreak\ \rho_{R}(y,c)}\Bigg{]}.

(78)

We now analyze Eq. (78) for both small and large rescaled temperatures $c$ . In the following, we use the value of chemical potential $\mu_{R}(c)$ which is obtained by numerical solving Eq. (78) with the constraint that the density is normalized to unity.

B.1.1 Small rescaled temperature $c\ll 1$

At zero temperature i.e., $c=0$ , all the hard rods arrange themselves leaving no gaps. In other words the center to center distance between the rods is $a$ , thereby making the density $\rho_{N}(x,0)=N/a$ where we recall that $N$ is number of hard rods. In the rescaled density variables this corresponds to a scaled density profile

\displaystyle\rho_{R}(y,0)=\begin{cases}\frac{1}{a}&\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ y\leq\left|\frac{a}{2}\right|\\ 0&\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ y>\left|\frac{a}{2}\right|.\end{cases}

(79)

We now study the effects of turning on a small temperature. More precisely we address how the zero temperature profile given in Eq. (79) gets smeared. Note that at

\displaystyle y=y_{c}=\left(\mu_{c}\delta\right)^{\frac{1}{\delta}}.

(80)

the square bracket in Eq. (78) changes sign. This in turn determines the following three distinct regions

(i)

Bulk region ( $|y|<y_{c}$ ): The density profile deviates from the value $1/a$ .
(ii)

Edge region (a zone where $|y-y_{c}|\lesssim O(c)$ ): The density profile deviates from a value $\rho_{R}(y_{c},c)=\rho^{*}_{R}$ which is the density at $y=y_{c}$ . At this value of $\rho^{*}_{R}$ the square bracket in Eq. (78) becomes zero.
(iii)

Tail region ( $|y|>y_{c}$ ): The density profile for finite temperature in this region deviates from its zero temperature value $\rho_{R}(y,0)=0$ .

We now compute the density profile at low temperatures of these three regions separately.

(i) Bulk region $|y|<y_{c}$ : In this region, we assume that the density takes the form

\displaystyle\rho_{R}(y,c)=\frac{1}{a}+\rho_{1}(y,c),

(81)

where $\rho_{1}(y,c)$ denotes the deviation about the zero temperature density. For sake of brevity, we henceforth omit the arguments of $\rho_{1}(y,c)$ . Using Eq. (81) in Eq. (78) we get

\displaystyle\mu_{c}-\frac{y^{\delta}}{\delta}

\displaystyle=-c\Bigg{[}\log{\left(\frac{-a^{2}\rho_{1}}{1+a\leavevmode\nobreak\ \rho_{1}}\right)}+\frac{1}{a\leavevmode\nobreak\ \rho_{1}}\Bigg{]}.

(82)

It turns out that a convenient perturbation parameter is the following

\displaystyle\nu(y)=\frac{c\delta}{\left(y_{c}^{\delta}-y^{\delta}\right)},

(83)

where we have used Eq. (80). Using Eq. (83) in Eq. (82) we obtain

\displaystyle-\frac{1}{\nu(y)}

\displaystyle=\log{\left(\frac{-a^{2}\rho_{1}}{1+a\leavevmode\nobreak\ \rho_{1}}\right)}+\frac{1}{a\leavevmode\nobreak\ \rho_{1}}.

(84)

To solve Eq. (84) we first perform a Taylor expansion

\displaystyle-\frac{1}{\nu(y)}

\displaystyle=\log{\left(-a^{2}\rho_{1}\right)}-a\leavevmode\nobreak\ \rho_{1}-\frac{a^{2}}{2}\rho_{1}^{2}+\frac{1}{a\leavevmode\nobreak\ \rho_{1}},

(85)

which can be again rearranged to give

\displaystyle a\leavevmode\nobreak\ \rho_{1}

\displaystyle=-\nu(y)\frac{1}{1+\nu(y)\left(\log{\left[-a^{2}\rho_{1}\right]}-a\leavevmode\nobreak\ \rho_{1}-\frac{a^{2}}{2}\rho_{1}^{2}\right)}.

(86)

We perform a Taylor series expansion [upto second order in $\nu(y)$ ] of the fraction on the right hand side of Eq. (86), since $\nu(y)\ll 1$ . This gives

	$\displaystyle\rho_{1}$	$\displaystyle\approx-\frac{\nu(y)}{a}\Bigg{[}1-\nu(y)\left(\log{\left[-a^{2}\rho_{1}\right]}-a\leavevmode\nobreak\ \rho_{1}-\frac{a^{2}}{2}\rho_{1}^{2}\right)$
		$\displaystyle+\nu(y)^{2}\left(\log{\left[-a^{2}\rho_{1}\right]}-a\leavevmode\nobreak\ \rho_{1}-\frac{a^{2}}{2}\rho_{1}^{2}\right)^{2}\Bigg{]}.$		(87)

Since the correction to the zero temperature density $\rho_{1}\ll 1$ and $\nu(y)\ll 1$ , we can invert Eq. (B.1.1) to express $\rho_{1}$ as a function of $\nu(y)$ order by order. This gives

$\displaystyle\rho_{1}^{(0)}$	$\displaystyle=-\frac{\nu(y)}{a},$	(88)
$\displaystyle\rho_{1}^{(1)}$	$\displaystyle=-\frac{\nu(y)}{a}+\frac{\nu(y)^{2}}{a}\log\Big{[}a\nu(y)\Big{]},$	(89)
$\displaystyle\rho_{1}^{(2)}$	$\displaystyle=-\frac{\nu(y)}{a}+\frac{\nu(y)^{2}}{a}\log\Big{[}a\nu(y)\Big{]}-\frac{\nu(y)^{3}}{a}\log\Big{[}a\nu(y)\Big{]}^{2},$	(90)

where recall that $\nu(y)$ is defined in Eq. (83). The superscript associated with $\rho_{1}$ in Eqs. (88)-(90) represents their respective orders. In Figs. 3a and 3d using Eqs. (88)-(90), we find a good agreement between the analytically obtained series solutions and the numerical solution of Eq. (78).

(ii) Edge region $|y-yc|\lesssim O(c)$ : Recall that this region is a zone defined by $|y-yc|\lesssim O(c)$ . Here $\nu(y)\gtrsim O(1)$ , and therefore, the above expressions Eqs. (88)- (90) fail. Hence, in this zone [ $|y_{c}-y|\lesssim O(c)$ ], we assume that the density takes the form

\displaystyle\rho_{R}(y,c)=\rho_{R}^{*}+\phi(y),

(91)

where $\rho_{R}^{*}$ is the value of the density at $y=y_{c}$ and the correction $\phi(y)\ll 1$ . The value of $\rho_{R}^{*}$ can be obtained by numerically solving Eq. (78) at $y=y_{c}$ which gives

\displaystyle\log{\left(\frac{1-a\rho_{R}^{*}}{\rho_{R}^{*}}\right)}-\frac{1}{1-a\rho_{R}^{*}}=0.

(92)

In this region, we define the perturbation parameter

\displaystyle b(y)=1-\exp\left[-\frac{y_{c}^{\delta}-y^{\delta}}{c\delta}\right]\ll 1.

(93)

Using Eq. (93) in Eq. (78) we get

	$\displaystyle\log{\left[1-b(y)\right]}$	$\displaystyle=\log{\Bigg{[}\frac{1-a\rho_{R}^{}-a\phi(y)}{\rho_{R}^{}+\phi(y)}\Bigg{]}}$
		$\displaystyle-\frac{1}{1-a\rho_{R}^{*}-a\phi(y)},$		(94)

which upon Taylor series expansion, upto third order in $\phi(y)$ , yields

$\displaystyle b(y)$	$\displaystyle\approx\phi(y)\left(\frac{1}{\rho_{R}^{}\left(1-a\rho_{R}^{}\right)^{2}}\right)$
	$\displaystyle-\phi(y)^{2}\left(\frac{-1+3a\rho_{R}^{}}{2\rho_{R}^{2}(1-a\rho_{R}^{*})^{3}}\right)$
	$\displaystyle-\phi(y)^{3}\left(\frac{-1+4a\rho_{R}^{}-6(a\rho_{R}^{})^{2}}{3\rho_{R}^{3}(1-a\rho_{R}^{})^{4}}\right).$	(95)

We can represent the correction to density $\phi(y)$ as a function of $b(y)$ by inverting Eq. (B.1.1) order by order which gives

$\displaystyle\phi^{(0)}(y)$	$\displaystyle=b(y)\rho_{R}^{}(1-a\rho_{R}^{})^{2},$	(96)
$\displaystyle\phi^{(1)}(y)$	$\displaystyle=b(y)\rho_{R}^{}(1-a\rho_{R}^{})^{2}$
	$\displaystyle+\frac{b(y)^{2}}{2}\rho_{R}^{}(1-a\rho_{R}^{})^{3}\big{(}-1+3a\rho_{R}^{*}\big{)},$	(97)
$\displaystyle\phi^{(2)}(y)$	$\displaystyle=b(y)\rho_{R}^{}(1-a\rho_{R}^{})^{2}$
	$\displaystyle+\frac{b(y)^{2}}{2}\rho_{R}^{}(1-a\rho_{R}^{})^{3}\big{(}-1+3a\rho_{R}^{*}\big{)}$
	$\displaystyle+\frac{b(y)^{3}}{6}\rho_{R}^{}(1-a\rho_{R}^{})^{4}\Big{(}1-10a\rho_{R}^{}+15(a\rho_{R}^{})^{2}\Big{)}.$	(98)

In Figs. 3b and 3e, using Eqs. (96)-(98), we compare the analytically obtained series solutions with the numerical solution of Eq. (78) and see reasonable agreement.

(iii) Tail region $|y|>y_{c}$ : In this region we assume that the density is very small and takes the form $\rho_{R}(y,c)=\rho_{1}$ with $\rho_{1}\ll 1$ . Using this assumption in Eq. (78) we get

\displaystyle\mu_{c}

\displaystyle=\frac{y^{\delta}}{\delta}-c\Bigg{[}\log{\left(\frac{1-a\leavevmode\nobreak\ \rho_{1}}{\rho_{1}}\right)}-\frac{1}{1-a\leavevmode\nobreak\ \rho_{1}}\Bigg{]}.

(99)

We introduce the perturbation parameter

\displaystyle\epsilon(y)=\exp\left(\frac{y_{c}^{\delta}-y^{\delta}}{c\delta}\right).

(100)

Since in this region $|y|>y_{c}$ and $c\ll 1$ , it implies $\epsilon(y)\ll 1$ . Using Eq. (100) in Eq. (78) we get

\displaystyle\log{\big{[}\epsilon(y)\big{]}}

\displaystyle=-\log{\left(\frac{1-a\leavevmode\nobreak\ \rho_{1}}{\rho_{1}}\right)}+\frac{1}{1-a\leavevmode\nobreak\ \rho_{1}}.

(101)

After some algebra and assuming $\rho_{1}\ll 1$ in Eq. (101), we obtain the transcendental equation

\displaystyle\rho_{1}

\displaystyle\approx\frac{\epsilon(y)}{e}\exp\Bigg{(}-2a\rho_{1}-\frac{(a\rho_{1})^{2}}{2}\Bigg{)},

(102)

where $e\approx 2.71828$ is the Euler’s number. We can now represent the density in terms of $\epsilon(y)/e$ as a series given by

$\displaystyle\rho_{1}^{(0)}$	$\displaystyle=\frac{\epsilon(y)}{e},$	(103)
$\displaystyle\rho_{1}^{(1)}$	$\displaystyle=\frac{\epsilon(y)}{e}-2a\left(\frac{\epsilon(y)}{e}\right)^{2},$	(104)
$\displaystyle\rho_{1}^{(2)}$	$\displaystyle=\frac{\epsilon(y)}{e}-2a\left(\frac{\epsilon(y)}{e}\right)^{2}+\frac{11}{2}a^{2}\left(\frac{\epsilon(y)}{e}\right)^{3}.$	(105)

In Figs. 3c and 3f, using Eqs. (103)-(105), we show a good agreement between the analytically obtained series solutions with the numerical solution of Eq. (78). Recall that the chemical potential $\mu_{R}(c)$ in the perturbation parameter $\epsilon(y)$ given in Eq. (100) is obtained from the numerical solution of Eq. (78) along with the constraint of unit normalized density.

B.1.2 Large rescaled temperatures: $c\gg 1$

When the temperature is high the particles explore a larger region in space thereby diluting the system as a consequence of which we get $\rho_{R}(y,c)\ll 1$ . We introduce a convenient perturbation parameter

\displaystyle\eta(y)=\exp\left(\frac{\mu_{c}\delta-y^{\delta}}{c\delta}\right),

(106)

where $\eta(y)\ll 1$ , since the chemical potential [see Fig. 4 in the main text], obtained by numerically solving Eq. (78), is negative and diverges for $c\gg 1$ .

We use the approximation $\rho_{R}(y,c)\ll 1$ in Eq. (78), and use a procedure mathematically similar to the low temperature tail region (Appendix. B.1.1), to obtain the expressions

$\displaystyle\rho_{R}^{(0)}(y,c)$	$\displaystyle=\frac{\eta(y)}{e},$	(107)
$\displaystyle\rho_{R}^{(1)}(y,c)$	$\displaystyle=\frac{\eta(y)}{e}-2a\left(\frac{\eta(y)}{e}\right)^{2},$	(108)
$\displaystyle\rho_{R}^{(2)}(y,c)$	$\displaystyle=\frac{\eta(y)}{e}-2a\left(\frac{\eta(y)}{e}\right)^{2}+\frac{11}{2}a^{2}\left(\frac{\eta(y)}{e}\right)^{3}.$	(109)

Note that the superscript in Eqs. (107)-(109) represents the order in $\eta(y)$ . In Fig. 5, we see a decent agreement of the analytically obtained series solutions with the numerical solution of Eq. (78).

B.2 Asymptotic densities for hyperbolic Calogero model

In the following subsections, we compute the asymptotic densities for the hyperbolic Calogero model at low and high rescaled temperature $c$ using a similar approach as described above for the HR model in Appendix B.1. Here we recall that the saddle point equation for the HC model is

\displaystyle\mu_{C}(c)=\frac{y^{\delta}}{\delta}+3\zeta(2)\rho_{C}(y,c)^{2}+c\log\rho_{C}(y,c).

(110)

As before we analyze Eq. (110) for small and large $c$ . In the following, we use the value of chemical potential $\mu_{C}(c)$ which is obtained by numerically solving Eq. (110) with the constraint that the density is normalized to unity.

B.2.1 Small rescaled temperatures $c\ll 1$

For $c\ll 1$ the densities $\rho_{C}(y,c)$ are assumed to be a small deviation from the zero temperature which is obtained by solving Eq. (110) for $c=0$ . The density profile is then given by

\displaystyle\rho_{C}(y,0)=\begin{cases}A_{\delta}\left(l^{\delta}-y^{\delta}\right)^{\frac{1}{2}}&\text{for}\quad|y|<l\\ 0&\text{for}\quad|y|>l,\end{cases}

(111)

where

\displaystyle A_{\delta}=\left(3\delta\zeta(2)\right)^{-\frac{1}{2}}

(112)

and the edge of the support of the density is given by

\displaystyle l=\left(\frac{\delta}{2A_{\delta}{\rm Beta}\left(\frac{1}{\delta},\frac{3}{2}\right)}\right)^{\frac{2}{2+\delta}}.

(113)

Here

\displaystyle{\rm Beta}(x,y)=\int_{0}^{1}\leavevmode\nobreak\ dr\leavevmode\nobreak\ r^{x-1}(1-r)^{y-1},

(114)

is the Beta function. The zero temperature chemical potential is given by [Eq. (51) of main text]

\displaystyle\mu_{C}(0)

\displaystyle=\left(\frac{\pi}{2}\right)^{\frac{\delta}{\delta+2}}\left(\frac{\delta^{-1/\delta}\Gamma\left(\frac{3}{2}+\frac{1}{\delta}\right)}{\Gamma\left(1+\frac{1}{\delta}\right)}\right)^{\frac{2\delta}{\delta+2}}.

(115)

Similar to the HR model [Appendix B.1.1], at low temperatures the value of chemical potential $\mu_{C}(c)$ [in Eq. (110)], determines the (i) bulk $|y|<y_{c}$ , (ii) edge $|y-y_{c}|\lesssim O(c)$ and the (iii) tail regions $|y|>y_{c}$ , where

\displaystyle y_{c}=\left(\mu_{C}(c)\delta\right)^{\frac{1}{\delta}}.

(116)

We compute the density profile at low temperatures separately for the three regions.
(i) Bulk region $|y|<y_{c}$ : In this region we assume that the density takes the form

\displaystyle\rho_{C}(y,c)=\rho_{C}(y,0)+\rho_{1}(y,c),

(117)

where $\rho_{1}(y,c)$ is the correction to the zero temperature density. We use the following notations in the rest of the calculations

	$\displaystyle\rho_{1}(y,c)\equiv\rho_{1},$	$\displaystyle\quad\quad\rho_{C}(y,0)\equiv\rho_{0}$
	$\displaystyle\mu_{C}(c)\equiv\mu_{c},$	$\displaystyle\quad\quad\mu_{C}(0)\equiv\mu_{0}.$		(118)

Using Eq. (B.2.1) in Eq. (110) gives

	$\displaystyle\mu_{c}$	$\displaystyle=\frac{y^{\delta}}{\delta}+c\log\rho_{0}+c\log\left[1+\frac{\rho_{1}}{\rho_{0}}\right]$
		$\displaystyle+3\zeta(2)\left(\rho_{0}+\rho_{1}\right)^{2}.$		(119)

Since the temperature is low ( $c\ll 1$ ), the correction to the zero temperature density $\rho_{1}\ll\rho_{0}$ . Furthermore we introduce the perturbation parameter in the bulk region

\displaystyle\nu(y)=\frac{\mu_{0}+c\log\rho_{0}-\mu_{c}}{c+6\zeta(2)\rho_{0}^{2}}\ll 1.

(120)

For $c\ll 1$ , it turns out that $\mu_{c}$ and $\mu_{0}$ are very close, which implies $\nu(y)\ll 1$ and therefore a suitable perturbative parameter. Using Eq. (120) in Eq. (B.2.1) and expanding Logarithmic term up to $(\rho_{1}/\rho_{0})^{3}$ gives

\displaystyle\frac{\rho_{1}}{\rho_{0}}

\displaystyle\approx-\nu(y)-\frac{1}{2}\left(\frac{\rho_{1}}{\rho_{0}}\right)^{2}-\frac{1}{3}\left(\frac{\rho_{1}}{\rho_{0}}\right)^{3}\frac{c}{c+6\zeta(2)\rho_{0}^{2}}.

(121)

We can solve Eq. (121) iteratively which gives

$\displaystyle\frac{\rho_{1}^{(0)}}{\rho_{0}}$	$\displaystyle=-\nu(y),$	(122)
$\displaystyle\frac{\rho_{1}^{(1)}}{\rho_{0}}$	$\displaystyle=-\nu(y)-\frac{\nu(y)^{2}}{2},$	(123)
$\displaystyle\frac{\rho_{1}^{(2)}}{\rho_{0}}$	$\displaystyle=-\nu(y)-\frac{\nu(y)^{2}}{2}-\nu(y)^{3}\left(\frac{1}{2}-\frac{1}{3}\frac{c}{c+6\zeta(2)\rho_{0}^{2}}\right),$	(124)

where the superscript represents their respective orders. In Fig. 8a and 8d using the Eqs. (122)-(124) we find a good agreement between the analytically obtained series solution and the numerical solution of Eq. (110).

(ii) Edge region $|y-y_{c}|\lesssim O(c)$ : Recall that this region is a zone defined by $|y-y_{c}|\lesssim O(c)$ . Here $\nu(y)$ defined in Eq. (120) is no longer a small parameter and therefore ill-suited to be a perturbation parameter. Hence as in the case of HR model we assume the density to takes the form

\displaystyle\rho_{C}(y,c)=\rho_{C}^{*}+\phi(y),

(125)

where $\rho_{C}^{*}$ is the value of the density at $y=y_{c}$ and the deviation $\phi(y)\ll\rho_{C}^{*}$ . The value of $\rho_{C}^{*}$ can be obtained by numerically solving Eq. (110) at $y=y_{c}$ where $y_{c}$ is given in Eq. (116). This gives

\displaystyle 3\zeta(2)\rho_{C}^{*2}+c\log\rho_{C}^{*}=0.

(126)

We now introduce a perturbative parameter

\displaystyle b(y)=\frac{y_{c}^{\delta}-y^{\delta}}{\delta\left(c+6\zeta(2)\rho_{C}^{*2}\right)}\ll 1,

(127)

since in this region $|y-y_{c}|\lesssim O(c)$ and $c\ll 1$ . Substituting Eq. (125) and using Eq. (127) in Eq. (110) and expanding, we get

\displaystyle\frac{\phi(y)}{\rho_{C}^{*}}

\displaystyle=b(y)-\frac{1}{2}\left(\frac{\phi(y)}{\rho_{C}^{*}}\right)^{2}-\frac{1}{3}\left(\frac{\phi(y)}{\rho_{C}^{*}}\right)^{3}\frac{c}{c+6\zeta(2)\rho^{*2}}.

(128)

Using a similar iterative approach as before we can represent the correction to the zero temperature density $\phi(y)$ as

$\displaystyle\frac{\phi^{(0)}(y)}{\rho_{C}^{*}}$	$\displaystyle=b(y),$	(129)
$\displaystyle\frac{\phi^{(1)}(y)}{\rho_{C}^{*}}$	$\displaystyle=b(y)-\frac{b(y)^{2}}{2},$	(130)
$\displaystyle\frac{\phi^{(2)}(y)}{\rho_{C}^{*}}$	$\displaystyle=b(y)-\frac{b(y)^{2}}{2}+b(y)^{3}\left(\frac{1}{2}-\frac{1}{3}\frac{c}{c+6\zeta(2)\rho_{0}^{2}}\right).$	(131)

In Figs. 8b. and 8e, using Eqs. (129)-(131), we find a good agreement between the analytically obtained series solution and the numerical solution of Eq. (110).

(iii) Tail region $|y|>y_{c}$ : Unlike the bulk and the edge regions, in the tail region we assume that the density is small and takes the form $\rho_{C}(y,c)=\rho_{1}$ where $\rho_{1}\ll 1$ . Using this assumption in Eq. (110) we obtain

\displaystyle\mu_{c}

\displaystyle=\frac{y^{\delta}}{\delta}+c\log\rho_{1}+3\zeta(2)\rho_{1}^{2}.

(132)

We now introduce a suitable perturbation parameter

\displaystyle\epsilon(y)=\exp\left(\frac{y_{c}^{\delta}-y^{\delta}}{c\delta}\right).

(133)

In the tail region, since $|y|>y_{c}$ and $c\ll 1$ , the perturbation parameter $\epsilon(y)\ll 1$ . Using Eq. (133) in Eq. (132) and rearranging the terms gives

\displaystyle\rho_{1}

\displaystyle=\epsilon(y)\exp\Bigg{(}-\frac{3\zeta(2)\rho_{1}^{2}}{c}\Bigg{)}.

(134)

We can now represent the density in terms of $\epsilon(y)$ by using the iterative scheme, similar to bulk and edge regions, on Eq. (134). This then gives

$\displaystyle\rho_{1}^{(0)}=\epsilon(y)$	,	(135)
$\displaystyle\rho_{1}^{(1)}=\epsilon(y)$	$\displaystyle\left[1-\frac{3\zeta(2)}{c}\epsilon(y)^{2}\right],$	(136)
$\displaystyle\rho_{1}^{(2)}=\epsilon(y)$	$\displaystyle\Bigg{[}1-\frac{3\zeta(2)}{c}\epsilon(y)^{2}+\frac{5}{2}\left(\frac{3\zeta(2)}{c}\right)^{2}\epsilon(y)^{4}\Bigg{]}.$	(137)

In Fig. 8c and 8f, using Eqs. (135)-(137), we show the analytically obtained asymptotic densities of the matches well with the numerical solution of Eq. (110).

B.2.2 Large rescaled temperature: $c\gg 1$

Similar to the HR model (Appendix B.1.2), when the temperature is high the particles are spread out over a larger region of space hence diluting the system as a consequence of which we get $\rho_{C}(y,c)\ll 1$ . We again introduce a convenient perturbation parameter

\displaystyle\eta(y)=\exp\left(\frac{\mu_{c}\delta-y^{\delta}}{\delta c}\right).

(138)

Since the chemical potential (see Fig. 7 in the main text), obtained by numerically solving Eq. (110), is negative and diverges for $c\gg 1$ , the perturbation parameter becomes small i.e., $\eta(y)\ll 1$ . Using Eq. (138) along with the low density approximation in Eq. (110) and following a procedure mathematically similar to low temperature tail region (Appendix B.2.1) we obtain the expressions

$\displaystyle\rho_{C}^{(0)}(y,c)=\eta(y)$	,	(139)
$\displaystyle\rho_{C}^{(1)}(y,c)=\eta(y)$	$\displaystyle\left[1-\frac{3\zeta(2)}{c}\eta(y)^{2}\right],$	(140)
$\displaystyle\rho_{C}^{(2)}(y,c)=\eta(y)$	$\displaystyle\Bigg{[}1-\frac{3\zeta(2)}{c}\eta(y)^{2}+\frac{5}{2}\left(\frac{3\zeta(2)}{c}\right)^{2}\eta(y)^{4}\Bigg{]}.$	(141)

In Fig. 9, using Eqs. (139)-(141), we see a good agreement of the analytically obtained series solutions with the numerical solution of Eq. (110).

Appendix C Nearest neighbor harmonic chain

In this section, we derive the analytical results of the density profiles for the nearest neighbor harmonic chain in a quadratic trap. Recall that, the equilibrium density in the Toda model shows a distinctly different behavior when compared with the HR and HC models as shown in Fig. 10 of the main text. We find that, unlike HR and HC models, the spatial spread of the equilibrium density profile is $N$ independent. We suspect that this peculiarity is rooted in the fact that the inter-particle repulsion is weak in the Toda model, thereby allowing the particle trajectories to cross. To elucidate the effects of weak inter-particle repulsion in the Toda model, we study the harmonic (Hessian) approximation of the Toda interaction $V_{T}(r_{i})$ [Eq. (56) of the main text], given by

\displaystyle V_{H}(r_{i})=\frac{k_{2}}{2}r_{i}^{2}-k_{3}r_{i},

(142)

where $r_{i}=x_{i+1}-x_{i}$ , $J>0$ (we set $J=1$ ), $k_{2}=1/d^{2}$ , and $k_{3}=1/d$ are the interaction strengths and $d$ determines the length scale of the interaction of the Toda model [Eq. (56) of the main text]. The subscript $H$ in Eq. (142) stands for the harmonic chain. Since the harmonic chain in quadratic trap is analytically tractable, it provides a transparent way for understanding the behavior of the Toda model in a suitable parameter regime.

In the following section, we find the equilibrium density profile of the harmonic oscillator in the quadratic trap which has the Hamiltonian

\displaystyle H=\frac{k_{1}}{2}\sum_{i=1}^{N}x_{i}^{2}+\frac{k_{2}}{2}\sum_{i=1}^{N-1}(x_{i+1}-x_{i})^{2}+k_{3}(x_{1}-x_{N}).

(143)

Here $k_{1}$ is the strength of the quadratic trap, $k_{2}$ is the spring constant and $k_{3}$ is the magnitude of the effective external force. The average thermal density is given by

\displaystyle\rho_{H}(x,\beta)=\frac{1}{N}\sum_{l=1}^{N}\langle\delta(x-x_{l})\rangle_{\beta}=\frac{1}{N}\sum_{l=1}^{N}P(x_{l}),

(144)

where the $\langle.\rangle_{\beta}$ is the thermal average in the canonical ensemble and $P(x_{l})$ is the marginal distribution of the position of the $l^{\rm th}$ particle. The task then is to find the marginal distribution, which is

	$\displaystyle P(x_{l}=x)=\frac{1}{Z_{H}}\int_{-\infty}^{\infty}dx_{1}..dx_{l}..dx_{N}$	$\displaystyle\leavevmode\nobreak\ \delta(x_{l}-x)\times$
		$\displaystyle{\rm exp}\left(-\beta H\big{(}\{x_{i}\}\big{)}\right),$		(145)

where the partition function is given by

\displaystyle Z_{H}=\int_{-\infty}^{\infty}dx_{1}..dx_{l}..dx_{N}\leavevmode\nobreak\ {\rm exp}\left(-\beta H\big{(}\{x_{i}\}\big{)}\right).

(146)

Using the variable transformation $y_{l}=\sqrt{\beta k_{2}}\leavevmode\nobreak\ x_{l}$ , Eq. (C) becomes

	$\displaystyle P\Bigg{(}x_{l}=\frac{y}{\sqrt{\beta k_{2}}}\Bigg{)}=\frac{\big{(}\beta k_{2}\big{)}^{\frac{1-N}{2}}}{Z_{H}}$	$\displaystyle\int_{-\infty}^{\infty}dy_{1}..dy_{l}..dy_{N}\leavevmode\nobreak\ \delta(y_{l}-y)\times$
		$\displaystyle{\rm exp}\left[-\beta H\Bigg{(}\Big{\{}\frac{y_{i}}{\sqrt{\beta k_{2}}}\Big{\}}\Bigg{)}\right],$		(147)

We note that the delta function in Eq. (C) has the following Fourier representation

\displaystyle\delta(y_{l}-y)=\int_{-\infty}^{\infty}\leavevmode\nobreak\ dq\exp\big{[}-iq(y_{l}-y)\big{]}.

(148)

Also note that the Hamiltonian in Eq. (143) can be recast in the rescaled variables ${\bf y^{T}}=[y_{1}..y_{l}..y_{N}]$ as

	$\displaystyle H=\frac{1}{\beta}\Big{[}\frac{1}{2}{\bf y^{T}Ay}+{\bf b^{T}y}\Big{]}$	$\displaystyle=\frac{1}{\beta}\Bigg{[}\frac{k}{2}\sum_{i=1}^{N}y_{i}^{2}+\frac{1}{2}\sum_{i=1}^{N-1}(y_{i+1}-y_{i})^{2}$
		$\displaystyle+\gamma(y_{1}-y_{N})+iqy_{l}\Bigg{]}.$		(149)

Here

\displaystyle k=\frac{k_{1}}{k_{2}},\quad\gamma=\sqrt{\frac{k_{3}^{2}\beta}{k_{2}}}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\quad b_{j}=\gamma(\delta_{1,j}-\delta_{N,j})+iq\delta_{jl}

(150)

is the $j^{\rm th}$ element of ${\bf b}$ . Since the interactions are nearest neighbour, ${\bf A}=[A_{ij}]$ is a $N\times N$ tridiagonal matrix given by

{\bf A}=\begin{bmatrix}\upsilon-1&-1&0&...&0&0\\ -1&\upsilon&-1&...&0&0\\ 0&-1&\upsilon&...&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&...&\upsilon&-1\\ 0&0&0&...&-1&\upsilon-1\\ \end{bmatrix},

(151)

with

\displaystyle\upsilon=k+2\leavevmode\nobreak\ .

(152)

Using Eq. (148) and Eq. (C) in Eq. (C) we get

	$\displaystyle P\left(x_{l}=\frac{y}{\sqrt{\beta k_{2}}}\right)$	$\displaystyle=\frac{\big{(}\beta k_{2}\big{)}^{\frac{1-N}{2}}}{Z_{H}}\int_{-\infty}^{\infty}\leavevmode\nobreak\ dq\leavevmode\nobreak\ \exp(iqy)\times$
		$\displaystyle\int_{-\infty}^{\infty}dy_{1}..dy_{l}..dy_{N}\leavevmode\nobreak\ {\rm exp}\left(-\frac{1}{2}{\bf y^{T}Ay}-{\bf b^{T}y}\right).$		(153)

Eq. (C) is a multivariate Gaussian integral which gives

	$\displaystyle P\left(x_{l}=\frac{y}{\sqrt{\beta k_{2}}}\right)=\frac{\big{(}\beta k_{2}\big{)}^{\frac{1-N}{2}}}{Z_{H}}$	$\displaystyle\frac{(2\pi)^{\frac{N}{2}}}{\sqrt{A_{N}}}\int_{-\infty}^{\infty}\leavevmode\nobreak\ dq\leavevmode\nobreak\ \exp(iqy)\times$
		$\displaystyle{\rm exp}\left(-\frac{1}{2}{\bf b^{T}A^{-1}b}\right),$		(154)

where $A_{N}$ is the determinant of matrix ${\bf A}$ of size $N\times N$ . It turns out the integral over $q$ in Eq. (C) is Gaussian which can be performed to obtain the normalized distribution of $l^{\rm th}$ particle

\displaystyle P\left(x_{l}=x\right)=\frac{1}{\sqrt{2\pi{\rm Var}(x_{l})}}{\rm exp}\left(-\frac{\Big{(}x-\langle x_{l}\rangle_{\beta}\Big{)}^{2}}{2\leavevmode\nobreak\ {\rm Var}(x_{l})}\right),

(155)

with the mean position given by

\displaystyle\langle x_{l}\rangle_{\beta}=\frac{\gamma}{\sqrt{\beta k_{2}}}(A^{-1}_{Nl}-A^{-1}_{1l}),

(156)

and the variance given by

\displaystyle{\rm Var}(x_{l})=\langle x_{l}^{2}\rangle_{\beta}-\langle x_{l}\rangle_{\beta}^{2}=\frac{A_{ll}^{-1}}{\beta k_{2}}.

(157)

We now have to compute the diagonal elements $A^{-1}_{ll}$ , along with the elements in the first and last row of ${\bf A^{-1}}$ i.e., $A^{-1}_{1l}$ and $A^{-1}_{Nl}$ . This can be explicitly computed based on the approach in Ref. [40] which gives

$\displaystyle A^{-1}_{1l}$	$\displaystyle=\frac{B_{N-l}}{A_{N}},$	(158)
$\displaystyle A^{-1}_{Nl}$	$\displaystyle=\frac{B_{l-1}}{A_{N}},$	(159)
$\displaystyle A^{-1}_{ll}$	$\displaystyle=\frac{B_{l-1}B_{N-l}}{A_{N}},$	(160)

where recall that $A_{N}$ is the determinant of the ${\bf A}$ [Eq. (151)] and $B_{N}$ is the determinant of an $N$ dimensional square matrix

{\bf B}=\begin{bmatrix}\upsilon&-1&0&...&0&0\\ -1&\upsilon&-1&...&0&0\\ 0&-1&\upsilon&...&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&...&\upsilon&-1\\ 0&0&0&...&-1&\upsilon-1\\ \end{bmatrix}.

(161)

We find that the determinants $A_{N}$ and $B_{N}$ are related through the recursion relations

	$\displaystyle A_{N}$	$\displaystyle=(\upsilon-1)B_{N-1}-B_{N-2},$		(162)
	$\displaystyle B_{N}$	$\displaystyle=(\upsilon-1)C_{N-1}-C_{N-2}.$		(163)

Here $C_{N}$ is the determinant of an $N\times N$ tridiagonal matrix given by

{\bf C}=\begin{bmatrix}\upsilon&-1&0&...&0&0\\ -1&\upsilon&-1&...&0&0\\ 0&-1&\upsilon&...&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&...&\upsilon&-1\\ 0&0&0&...&-1&\upsilon\\ \end{bmatrix},

(164)

where recall that $\upsilon=k+2$ . The determinant of the $N$ -dimensional ${\bf C}$ matrix has been computed in Ref. [40] and is given by

\displaystyle C_{N}

\displaystyle=\frac{\chi^{-N}}{\chi^{2}-1}\left(\chi^{2N+2}-1\right),

(165)

where

\displaystyle\chi=\frac{\upsilon+\sqrt{\upsilon^{2}-4}}{2}.

(166)

Using the recursion relations Eqs. (162) and (163) with the determinant of $C$ matrix Eq. (165) we get

	$\displaystyle B_{N}$	$\displaystyle=\frac{\chi^{-N}}{\chi+1}\left(\chi^{2N+1}+1\right),$		(167)
	$\displaystyle A_{N}$	$\displaystyle=\frac{\chi^{-N}(\chi-1)}{\chi+1}\left(\chi^{2N}-1\right).$		(168)

Using these expressions, Eqs. (167) and (168), in Eqs. (156)-(160) we get

	$\displaystyle\langle x_{l}\rangle_{\beta}=$	$\displaystyle\frac{\gamma}{\sqrt{\beta k_{2}}}\frac{\left(\chi^{N+l}+\chi^{N-l+1}-\chi^{2N-l+1}-\chi^{l}\right)}{(\chi-1)(\chi^{2N}-1)},$		(169)
	$\displaystyle{\rm Var}(x_{l})$	$\displaystyle=\frac{\chi}{\beta k_{2}}\frac{\left(1+\chi^{2N}+\chi^{2N+1-2l}+\chi^{2l-1}\right)}{\left(\chi^{2}-1\right)\left(\chi^{2N}-1\right)}.$		(170)

Using the expression of the marginal distribution $P(x_{l})$ [Eq. (155)] in the expression for average density profile Eq. (144), we get the expression, for the case of harmonic chain in the quadratic trap

\displaystyle\rho_{H}\left(x,\beta\right)

\displaystyle=\frac{1}{N}\sum_{l=1}^{N}\frac{1}{\sqrt{2\pi{\rm Var}(x_{l})}}{\rm exp}\left(-\frac{\Big{(}x-\langle x_{l}\rangle_{\beta}\Big{)}^{2}}{2\leavevmode\nobreak\ {\rm Var}(x_{l})}\right).

(171)

Eq. (171) is the exact result for the harmonic chain in quadratic trap. We now compare Eq. (171) with the Monte-Carlo (MC) densities of the Toda model in quadratic trap in appropriate parameter regimes.

In Fig. 12, we show a comparison of the the MC density profiles of the Toda model, for $N=32,64,128$ , with the analytically obtained densities of the harmonic chain model [see Eq. (171)], for two values of the Toda interaction length scale $d=1.0,10.0$ . We find a good agreement between the densities of both the models which suggests the $N$ -independence of the density profiles. This $N$ -independence can be understood in the case of harmonic chain in quadratic trap as follows. For $d\gg 1$ , Eq. (166) can be approximated as

\displaystyle\chi=\frac{\upsilon+\sqrt{\upsilon^{2}-4}}{2}\approx d^{2},

(172)

where $\upsilon=2+k=2+d^{2}$ . Therefore, the mean [Eq. (156)] and variance [Eq. (157)] of the $l^{\rm th}$ particle position are given by

	$\displaystyle\langle x_{N}\rangle_{\beta}=-\langle x_{1}\rangle_{\beta}\approx\frac{1}{d}$	$\displaystyle\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \langle x_{l}\rangle_{\beta}\approx 0\leavevmode\nobreak\ \forall\leavevmode\nobreak\ l\in\{2..,N-1\}$		(173)
	$\displaystyle{\rm Var}(x_{l})$	$\displaystyle\approx\frac{1}{\beta}.$		(174)

Note that the mean and the variance of the position of the $l^{\rm th}$ particle given in Eqs. (173) and (174) respectively are independent of $N$ as can be seen from Figs. 14b and 14d. This proves that the density has no $N$ dependence for large $d$ and therefore Eq. (171) further simplifies to

\displaystyle\rho_{H}(x,\beta)\approx\sqrt{\frac{\beta}{2\pi}}\exp\Bigg{(}-\beta\frac{x^{2}}{2}\Bigg{)}.

(175)

These $N$ -independent density profiles are in accordance with our observations for both Toda and the harmonic chain models [see Figs. 10c, d and 12].

For small values of $d$ the system behaves as a very stiff harmonic chain since the spring constant $k_{2}=1/d^{2}$ is very large. For these values we expect that the thermodynamic limit can only be attained for extremely large values of $N$ . This can be understood more clearly when we consider the asymptotic behavior for the small values of $d$ . For small $d\ll 1$ , Eq. (166) becomes

\displaystyle\chi=\frac{\upsilon+\sqrt{\upsilon^{2}-4}}{2}\approx 1+d,

(176)

where $\upsilon=2+k=2+d^{2}$ . Using Eq. (176) in the expressions Eqs. (169) and (170), we obtain the mean and the variance of the $l^{\rm th}$ particle position

	$\displaystyle\langle x_{l}\rangle_{\beta}$	$\displaystyle=d\left(\frac{1}{(1+d)^{N-l+1}}-\frac{1}{(1+d)^{l}}\right),$		(177)
	$\displaystyle{\rm Var}(x_{l})$	$\displaystyle=\frac{d^{2}}{2\beta}\left(\frac{1}{1+d}+\frac{1}{(1+d)^{2l}}+\frac{1}{(1+d)^{2(N-l+1)}}\right),$		(178)

for $Nd\gg 1$ . In Fig. 14a and 14c, we show that the mean and the variance, computed using Eqs. (169) and Eq. (170) respectively, converges very slowly with increasing $N$ to zero and $d/(2\beta)$ respectively. Therefore, we expect the density profiles in Eq. (171) to slowly converge to an $N$ -independent profile for very large $N$ . The density profile in the large $N$ limit is then given by

\displaystyle\rho_{H}(x,\beta)\approx\sqrt{\frac{\beta}{d\pi}}\exp\Bigg{(}-\beta\frac{x^{2}}{d}\Bigg{)},

(179)

which is the curve corresponding to $N\to\infty$ in Fig. 13. This slow convergence is similar to what we find for the Toda model for $d=0.1$ as shown in Figs. 10b.

References

Babelon et al. [2003] O. Babelon, D. Bernard, and M. Talon, Introduction to classical integrable systems (Cambridge University Press, 2003).
Torrielli [2016] A. Torrielli, Classical integrability, J. Phys. A: Mathematical and Theoretical 49, 323001 (2016).
Cao [1995] C. Cao, Classical integrable systems, in Soliton Theory and Its Applications (Springer, 1995) pp. 152–191.
Bastianello et al. [2021] A. Bastianello, A. De Luca, and R. Vasseur, Hydrodynamics of weak integrability breaking, J. Stat. Mech.: Theory and Experiment 2021, 114003 (2021).
Kinoshita et al. [2006] T. Kinoshita, T. Wenger, and D. S. Weiss, A quantum newton’s cradle, Nature 440, 900 (2006).
Tang et al. [2018] Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan, and B. L. Lev, Thermalization near integrability in a dipolar quantum newton’s cradle, Phys. Rev. X 8, 021030 (2018).
Bouchoule and Dubail [2022] I. Bouchoule and J. Dubail, Generalized hydrodynamics in the one-dimensional bose gas: theory and experiments, J. Stat. Mech.: Theory and Experiment 2022, 014003 (2022).
Di Cintio et al. [2018] P. Di Cintio, S. Iubini, S. Lepri, and R. Livi, Transport in perturbed classical integrable systems: The pinned toda chain, Chaos, Solitons & Fractals 117, 249 (2018).
Lam and Kurchan [2014] K.-D. N. T. Lam and J. Kurchan, Stochastic perturbation of integrable systems: a window to weakly chaotic systems, J. Stat. Phys. 156, 619 (2014).
Rajabpour and Sotiriadis [2014] M. Rajabpour and S. Sotiriadis, Quantum quench of the trap frequency in the harmonic calogero model, Phys. Rev. A 89, 033620 (2014).
Cao et al. [2018] X. Cao, V. B. Bulchandani, and J. E. Moore, Incomplete thermalization from trap-induced integrability breaking: Lessons from classical hard rods, Phys. Rev. Lett. 120, 164101 (2018).
Dhar et al. [2019] A. Dhar, A. Kundu, J. L. Lebowitz, and J. A. Scaramazza, Transport properties of the classical toda chain: effect of a pinning potential, J. Stat. Phys. 175, 1298 (2019).
Caux et al. [2019] J.-S. Caux, B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, Hydrodynamics of the interacting bose gas in the quantum newton cradle setup, SciPost Physics 6, 070 (2019).
Durnin et al. [2021] J. Durnin, M. Bhaseen, and B. Doyon, Nonequilibrium dynamics and weakly broken integrability, Phys. Rev. Lett. 127, 130601 (2021).
Bulchandani et al. [2021] V. B. Bulchandani, M. Kulkarni, J. E. Moore, and X. Cao, Quasiparticle kinetic theory for calogero models, J. Phys. A: Mathematical and Theoretical 54, 474001 (2021).
Inguscio et al. [2008] M. Inguscio, W. Ketterle, and C. Salomon, Ultra-cold Fermi gases, Vol. 164 (IOS press, 2008).
Dalfovo et al. [1999] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of bose-einstein condensation in trapped gases, Reviews of Modern Physics 71, 463 (1999).
Joseph et al. [2011] J. Joseph, J. E. Thomas, M. Kulkarni, and A. G. Abanov, Observation of shock waves in a strongly interacting fermi gas, Phys. Rev. Lett. 106, 150401 (2011).
Tonks [1936] L. Tonks, The complete equation of state of one, two and three-dimensional gases of hard elastic spheres, Phys. Rev. 50, 955 (1936).
Olshanetsky and Perelomov [1981] M. A. Olshanetsky and A. M. Perelomov, Classical integrable finite-dimensional systems related to lie algebras, Physics Reports 71, 313 (1981).
Perelomov et al. [1990] A. M. Perelomov et al., Integrable Systems of Classical Mechanics and Lie Algebras: Volume I (Springer, 1990).
Polychronakos [1992] A. P. Polychronakos, New integrable systems from unitary matrix models, Phys. Lett. B 277, 102 (1992).
Fu et al. [2019] W. Fu, Y. Zhang, and H. Zhao, Universal law of thermalization for one-dimensional perturbed toda lattices, New Journal of Physics 21, 043009 (2019).
Gurappa and Panigrahi [1999] N. Gurappa and P. K. Panigrahi, Equivalence of the calogero-sutherland model to free harmonic oscillators, Phys. Rev. B 59, R2490 (1999).
Kulkarni and Polychronakos [2017] M. Kulkarni and A. Polychronakos, Emergence of the calogero family of models in external potentials: duality, solitons and hydrodynamics, J. Phys. A: Mathematical and Theoretical 50, 455202 (2017).
Gon and Kulkarni [2019] A. K. Gon and M. Kulkarni, Duality in a hyperbolic interaction model integrable even in a strong confinement: multi-soliton solutions and field theory, J. Phys. A: Mathematical and Theoretical 52, 415201 (2019).
Dean and Majumdar [2006] D. S. Dean and S. N. Majumdar, Large deviations of extreme eigenvalues of random matrices, Phys. Rev. Lett. 97, 160201 (2006).
Majumdar and Schehr [2014] S. N. Majumdar and G. Schehr, Top eigenvalue of a random matrix: large deviations and third order phase transition, J. Stat. Mech.: Theory and Experiment 2014, P01012 (2014).
Agarwal et al. [2019] S. Agarwal, A. Dhar, M. Kulkarni, A. Kundu, S. N. Majumdar, D. Mukamel, and G. Schehr, Harmonically confined particles with long-range repulsive interactions, Phys. Rev. Lett. 123, 100603 (2019).
Kumar et al. [2020] A. Kumar, M. Kulkarni, and A. Kundu, Particles confined in arbitrary potentials with a class of finite-range repulsive interactions, Phys. Rev. E 102, 032128 (2020).
Kethepalli et al. [2021] J. Kethepalli, M. Kulkarni, A. Kundu, S. N. Majumdar, D. Mukamel, and G. Schehr, Harmonically confined long-ranged interacting gas in the presence of a hard wall, J. Stat. Mech.: Theory and Experiment 2021, 103209 (2021).
Kethepalli et al. [2022] J. Kethepalli, M. Kulkarni, A. Kundu, S. N. Majumdar, D. Mukamel, and G. Schehr, Edge fluctuations and third-order phase transition in harmonically confined long-range systems, J. Stat. Mech.: Theory and Experiment 2022, 033203 (2022).
Santra et al. [2022] S. Santra, J. Kethepalli, S. Agarwal, A. Dhar, M. Kulkarni, and A. Kundu, Gap statistics for confined particles with power-law interactions, Phys. Rev. Lett. 128, 170603 (2022).
Mussardo [2010] G. Mussardo, Statistical field theory: an introduction to exactly solved models in statistical physics (Oxford University Press, 2010).
Cunden et al. [2018] F. D. Cunden, P. Facchi, M. Ligabò, and P. Vivo, Universality of the weak pushed-to-pulled transition in systems with repulsive interactions, J. Phys. A: Mathematical and Theoretical 51, 35LT01 (2018).
Hardin et al. [2018] D. P. Hardin, T. Leblé, E. B. Saff, and S. Serfaty, Large deviation principles for hypersingular riesz gases, Constructive Approximation 48, 61 (2018).
Allez et al. [2012] R. Allez, J.-P. Bouchaud, and A. Guionnet, Invariant beta ensembles and the gauss-wigner crossover, Phys. Rev. Lett. 109, 094102 (2012).
Percus [1976] J. Percus, Equilibrium state of a classical fluid of hard rods in an external field, J. Stat. Phys. 15, 505 (1976).
Stone et al. [2008] M. Stone, I. Anduaga, and L. Xing, The classical hydrodynamics of the calogero–sutherland model, J. Phys. A: Mathematical and Theoretical 41, 275401 (2008).
Hu and O’Connell [1996] G. Hu and R. F. O’Connell, Analytical inversion of symmetric tridiagonal matrices, J. Phys. A: Mathematical and General 29, 1511 (1996).

$\displaystyle b(y)$	$\displaystyle\approx\phi(y)\left(\frac{1}{\rho_{R}^{}\left(1-a\rho_{R}^{}\right)^{2}}\right)$
	$\displaystyle-\phi(y)^{2}\left(\frac{-1+3a\rho_{R}^{}}{2\rho_{R}^{2}(1-a\rho_{R}^{*})^{3}}\right)$
	$\displaystyle-\phi(y)^{3}\left(\frac{-1+4a\rho_{R}^{}-6(a\rho_{R}^{})^{2}}{3\rho_{R}^{3}(1-a\rho_{R}^{})^{4}}\right).$	(95)

Finite temperature equilibrium density profiles of integrable systems in confining potentials

Abstract

I Introduction

II Models and definitions

III Field theory formalism

IV Results from field theory and comparison with Monte-Carlo simulations

IV.1 Hard rods model

IV.2 Hyperbolic Calogero model

IV.3 Integrable models with particle-crossing

V Conclusion

Acknowledgements

Appendix A Derivation of the free energy

A.1 Free energy for hard rods

A.2 Free energy for hyperbolic Calogero model

Appendix B Analytical forms of density profiles for hard rods and hyperbolic Calogero model at low and high rescaled temperatures cc

B.1 Analytical forms of the density profiles for hard rods

B.1.1 Small rescaled temperature c≪1c\ll 1

B.1.2 Large rescaled temperatures: c≫1c\gg 1

B.2 Asymptotic densities for hyperbolic Calogero model

B.2.1 Small rescaled temperatures c≪1c\ll 1

B.2.2 Large rescaled temperature: c≫1c\gg 1

Appendix C Nearest neighbor harmonic chain

References

Appendix B Analytical forms of density profiles for hard rods and hyperbolic Calogero model at low and high rescaled temperatures $c$

B.1.1 Small rescaled temperature $c\ll 1$

B.1.2 Large rescaled temperatures: $c\gg 1$

B.2.1 Small rescaled temperatures $c\ll 1$

B.2.2 Large rescaled temperature: $c\gg 1$