# Supplementary material for ‘Gap statistics for confined particles with
power-law interactions’

S. Santra<EMAIL_ADDRESS>International Centre for Theoretical
Sciences, Tata Institute of Fundamental Research, Bengaluru – 560089, India
J. Kethepalli<EMAIL_ADDRESS>International Centre for
Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru –
560089, India  S. Agarwal<EMAIL_ADDRESS>Department of Physics,
University of Colorado, Boulder, Colorado 80309, USA  A. Dhar
<EMAIL_ADDRESS>International Centre for Theoretical Sciences, Tata
Institute of Fundamental Research, Bengaluru – 560089, India  M. Kulkarni
<EMAIL_ADDRESS>International Centre for Theoretical Sciences, Tata
Institute of Fundamental Research, Bengaluru – 560089, India  A. Kundu
<EMAIL_ADDRESS>International Centre for Theoretical Sciences, Tata
Institute of Fundamental Research, Bengaluru – 560089, India

###### Contents

  1. I Model: Riesz gas and density profile
  2. II Variance of the bulk gap
    1. II.1 Low temperature dynamics
    2. II.2 Details of Monte-Carlo (MC) simulations
    3. II.3 Microscopic Hessian
    4. II.4 Quadratic field theory
    5. II.5 System size scaling of the matrix elements $M$ and $H$
    6. II.6 Error estimation in the exponent $b_{k}$
  3. III Interesting aspects of the gap statistics for $-1\leq k\leq 0$
  4. IV Calogero-Moser system

This supplementary is organized as follows. In Sec. I, we give the details of
the model and recap some relevant previous results. In Sec. II we give details
of the computation of the variance of the bulk gap using Monte-Carlo
simulation (MC), microscopic Hessian (MH) and quadratic field theory (FT). In
Sec. III some interesting features of gap statistics for $-1\leq k\leq 0$ are
discussed. In Sec. IV, we study the gap distribution for the Calogero-Moser
(CM) model ($k=2$) and discuss reasons for disagreement with previous studies
Bogomolny  _et al._ (2009).

##  I Model: Riesz gas and density profile

In the main text we considered a classical system of $N$ particles on a line
confined by a harmonic potential and interacting with each other via a
repulsive power law interaction of the form $\text{sgn}(k)|r|^{-k}$, where $r$
is the distance between two particles. The potential energy of the system is
given by ($\forall k>-2$) Riesz (1938)

$\displaystyle
E\left(\\{x_{i}\\}\right)=\sum_{i=1}^{N}\frac{x^{2}_{i}}{2}+\frac{J~{}\text{sgn}(k)}{2}\sum_{i\neq
j}\frac{1}{\mid x_{i}-x_{j}\mid^{k}},$ (S1)

where $J>0$ and $\text{sgn}(k)$ ensures a repulsive interaction. Using large
$N$ field theory, the average thermal density is calculated recently Agarwal
_et al._ (2019) $\forall$ $k>-2$. The average density, $\rho^{\rm(eq)}_{N}$,
for large $N$ is independent of the inverse temperature, $\beta$ and is
described by the following scaling form

$\rho^{\rm(eq)}_{N}(x)\equiv\frac{1}{N}\sum_{j=1}^{N}\langle\delta(x-x_{j})\rangle=\frac{1}{l_{k}N^{\alpha_{k}}}F_{k}\big{(}\frac{x}{l_{k}N^{\alpha_{k}}}\big{)},$
(S2)

where $\langle...\rangle$ denotes an equilibrium average and the exponent
$\alpha_{k}$ is given by

$\displaystyle\alpha_{k}=\begin{cases}\frac{k}{k+2}~{}\rm{for}~{}k>1\\\
\frac{1}{k+2}~{}\rm{for}~{}-2<k<1.\end{cases}$ (S3)

The scaling function $F_{k}(y)$, supported over $y\in[-1/2,1/2]$, is given
explicitly by

$F_{k}(y)=\frac{1}{B\big{(}\gamma_{k}+1,\gamma_{k}+1\big{)}}\big{(}\frac{1}{4}-y^{2}\big{)}^{\gamma_{k}},$
(S4)

where the exponent $\gamma_{k}$ is given by

$\displaystyle\gamma_{k}=\begin{cases}\frac{1}{k}~{}~{}~{}~{}\rm{for}~{}k>1\\\
\frac{k+1}{2}~{}\rm{for}~{}-2<k<1,\end{cases}$ (S5)

and the system size independent length scale is explicitly

$\displaystyle
l_{k}=\begin{cases}\Big{(}\frac{\big{(}2J\zeta(k)(k+1)\big{)}^{1/k}}{B(1+1/k,1+1/k)}\Big{)}^{\frac{k}{k+2}}~{}~{}~{}~{}~{}~{}~{}~{}~{}\rm{for}~{}k>1\\\
\Big{(}\frac{J|k|\pi(k+1)}{\sin\big{[}\frac{\pi}{2}(k+1)\big{]}B(\frac{k+3}{2},\frac{k+3}{2})}\Big{)}^{\frac{1}{k+2}}~{}\rm{for}~{}-2<k<1,\end{cases}$
(S6)

where $B$ is beta function. We use these densities to study the statistical
properties of interparticle bulk gaps in detail.

##  II Variance of the bulk gap

In this section we discuss the calculation of the variance of bulk gap
$\Delta_{i}=x_{i+1}-x_{i}$ ($1\ll i\ll N-1$), where $\\{x_{i}\\}$ are the
positions of the particles. As discussed in the main text, the mean of the
bulk gap scales with system size as $\langle\Delta_{i}\rangle\sim
N^{\alpha_{k}}/N=N^{-a_{k}}$ where $a_{k}=1-\alpha_{k}$, i.e.,

$\displaystyle a_{k}=\begin{cases}\frac{2}{k+2}&\text{for }k>1\\\
\frac{k+1}{k+2}&\text{for}-2<k<1.\end{cases}$ (S7)

Our conjectured form for the variance of bulk gap is
$\sigma^{2}_{\Delta_{i}}\sim N^{-b_{k}}$, with

$\displaystyle b_{k}=\begin{cases}2~{}~{}&{\rm{for}}~{}~{}k>1\\\ 1+k&{\rm
for}~{}0<k<1\\\ 2(k+1)/(k+2)&{\rm for}~{}-1<k<0\\\ 1+k~{}~{}&{\rm
for}~{}~{}-2<k<-1.\end{cases}$ (S8)

Here we present the details of the microscopic Hessian (MH) and quadratic
field theory (FT) which are used to compute the variance and provide support
for the above conjectured form.

###  II.1 Low temperature dynamics

Figure S1: Spacetime trajectories of particles evolving under overdamped
Langevin equations and starting from their equilibrium positions
$\\{\bar{y}_{i}\\}$. The x-axis represents the scaled position
$z=x/(l_{k}N^{\alpha_{k}})$ and we show results for different values of $k$.
We took $\gamma=1,\beta=16$ and $N=64$.

To get some intuition about the microscopic fluctuations of particle positions
in our many-body system, it is instructive to study the overdamped Langevin
dynamics of the particles at low temperatures
($T=\beta^{-1}<N^{2\alpha_{k}}$). We thus consider the equations of motion

$\dot{x_{i}}(t)=-{\partial_{x_{i}}{E}}/{\gamma}+\sqrt{(2/\gamma\beta)}\eta_{i}(t),$
(S9)

where $\eta_{i}(t)$ are Gaussian white noise terms with zero mean and unit
variance. The initial condition are taken to be $x_{i}(0)=\bar{y}_{i}$, where
$\bar{y}_{i}$ are the equilibrium positions. This is computed from the
equilibrium density profile using

$i=N\int_{-l_{k}}^{\bar{w}_{i}}\rho^{\rm(eq)}_{N}(x)dx,$ (S10)

and then taking the symmetrised form

$\bar{y}_{i}=\frac{\bar{w}_{i}-\bar{w}_{N-i+1}}{2}.$ (S11)

From the spacetime trajectories shown in Fig. S1 we see that for all $k>-2$,
the particles fluctuate around their equilibrium positions. We find that the
fluctuations about the equilibrium configurations are large for $-1<k<0$
compared to other values of $k$ — this is consistent with our findings on the
gap flutuations relative to mean (see Sec. III).

###  II.2 Details of Monte-Carlo (MC) simulations

We compute the mean and variance of the bulk gap using Monte-Carlo
simulations. For a given processor, we disregard about initial $10^{6}$ MC
cycles. By MC cycle, we mean $O(N)$ Metropolis steps. We then collect data of
about $10^{6}$ samples. This process is simultaneously executed in around
$N/2$ processors thereby making number of samples to be around $10^{8}$.
Needless to mention, the exact numbers depend on whether we achieve
satisfactory convergence of our results for a given $N$ and $k$. For all the
spacing distributions (which contain a total of $N-1$ gaps ), we have used
$10^{6}$ MC samples. Both $J_{0}$ and $J$ are taken to be unity for log-gas
($k=0^{+}$) and for any other values of $k>-2$ respectively.

###  II.3 Microscopic Hessian

At equilibrium the particles are sampled from the Gibbs Boltzmann distribution

$P_{G}(\\{x_{i}\\})=\frac{e^{-\beta E\\{x_{i}\\}}}{Z},$ (S12)

where the partition function $Z=\int\prod_{i}dx_{i}e^{-\beta E\\{x_{i}\\}}$.
Then the joint distribution of all gaps is given by

$\mathcal{P}(\\{\Delta_{i}\\})=\int\left(\prod_{i=1}^{N-1}dx_{i}\delta\left(\Delta_{i}-\left(x_{i+1}-x_{i}\right)\right)\right)e^{-\beta
E(\\{x_{i}\\})}.$ (S13)

Using this distribution, in principle one can calculate mean and the variance
of the $i^{\rm th}$ gap. However the integrals over microscopic positions are
difficult to compute, so we only solve them approximately. At zero
temperature, the system will be in the ground state characterised by the
configuration of positions $y_{i}$ and corresponding gaps $\Delta_{i}^{\rm
GS}=y_{i+1}-y_{i}$ . At low temperatures we expect that the Hessian of the
microscopic Hamiltonian Eq. (S1) about the ground state would approximately
capture the behavior of the flucutations of the gap. The positions $y_{i}$ are
obtained by minimising the energy using the Broyden-Fletcher-Goldfarb-Shanno
(BFGS) algorithm Fletcher (1970); Broyden (1970). Under the Hessian
approximation the joint distribution of fluctuation of gaps
$\delta\Delta_{i}=\Delta_{i}-\Delta^{\rm GS}_{i}$

$\mathcal{P}_{MH}(\\{\Delta_{i}\\})\sim
e^{-\beta\left[\frac{1}{2}\sum_{i,j=1}^{N}H_{ij}\delta\Delta_{i}\delta\Delta_{j}\right]},$
(S14)

where the Hessian of the system about the equilibrium configuration is
$H_{ij}=\Big{[}\frac{\partial^{2}{E}}{\partial{\Delta_{i}}\partial{\Delta_{j}}}\Big{]}_{\rm
GS}$. The matrix $H_{ij}$ is obtained as follows. The Hessian $\tilde{H}_{ij}$
of the Hamiltonian Eq. (S1) in terms of the fluctuation of positions around
their equilibrium values is given by

$\begin{split}\tilde{H}_{ij}&=\Big{[}\frac{\partial^{2}{E(\\{x\\})}}{\partial{x_{i}}\partial{x_{j}}}\Big{]}_{\rm
GS}\\\ &=\delta_{ij}\Big{[}1+\sum_{n\neq
i}^{N}\frac{J\text{sgn}(k)k(k+1)}{\big{(}y_{i}-y_{n}\big{)}^{k+2}}\Big{]}\\\
&\quad\quad-(1-\delta_{ij})\frac{J\text{sgn}(k)k(k+1)}{\big{(}y_{i}-y_{j}\big{)}^{k+2}}.\end{split}$
(S15)

To get the Hessian $H_{ij}$ we change position to the gap variables,
$\Delta_{i}=x_{i+1}-x_{i}$ for $i=1,\ldots,N-1$ and the centre of mass
coordinate which we denote by $\Delta_{N}\equiv\sum_{i=1}^{N}x_{i}/N$ by the
transformation

$\displaystyle x_{j}$
$\displaystyle=\Delta_{N}-\sum_{i=j}^{N-1}\Delta_{i}+\sum_{i=1}^{N-1}\frac{i}{N}\Delta_{i}$
(S16) $\displaystyle=\sum_{j=1}^{N}A_{ij}\Delta_{j},$ (S17)

where $A$ is a matrix of dimension $N\times N$ with matrix elements

$\displaystyle A_{ij}=\begin{cases}\frac{j}{N}&\text{for }i>j\\\
(\frac{j}{N}-1)&\text{for }i\leq j(\neq N)\\\ 1&\text{for
}j=N,i=1,...,N.\end{cases}$ (S18)

The quadratic Hamiltonian in terms of new variables becomes
$\frac{1}{2}\sum_{i,j=1}^{N}H_{ij}\delta\Delta_{i}\delta\Delta_{j}$, where
$H=A^{T}\tilde{H}A$. Using Eq. (S14) we compute the variance of bulk gap from
the relation $\sigma^{2}_{i}=\langle\delta\Delta^{2}_{i}\rangle={(\beta
H)^{-1}_{ii}},$ where $H^{-1}_{ii}$ is found numerically.

###  II.4 Quadratic field theory

As discussed in Ref. Agarwal  _et al._ (2019) the Reisz gas for large-$N$ can
be described by a free energy functional
$\Sigma[\rho_{N}]=\mathcal{E}[\rho_{N}]-\beta^{-1}S[\rho_{N}]$ corresponding
to a macroscopic density profile $\rho_{N}(x)$ where $\mathcal{E}[\rho_{N}]$
is the energy and

$\displaystyle S[\rho_{N}]=-N\int_{-\infty}^{\infty}dx\rho_{N}\log(\rho_{N})$
(S19)

is the entropy functional. The form of the energy functional depends on $k$,
being local for $k\geq 1$ and non-local for $-2<k<1$ Agarwal  _et al._ (2019).
The energy functional,$\mathcal{E}[\rho_{N}(x)]$ (Agarwal  _et al._ (2019)) is
given by

$\displaystyle\mathcal{E}[\rho_{N}(x)]\approx$
$\displaystyle\frac{N}{2}\int_{-\infty}^{\infty}x^{2}\rho_{N}(x)dx-\mu\left(\int_{-\infty}^{\infty}\rho_{N}(x)dx-1\right)$
$\displaystyle+J\zeta(k)N^{k+1}\int_{-\infty}^{\infty}[\rho_{N}(x)]^{k+1}dx,$
(S20)

for $k>1$ whereas for $-2<k<1$, it is given by

$\displaystyle\mathcal{E}[\rho_{N}(x)]\approx\frac{N}{2}\int_{-\infty}^{\infty}x^{2}\rho_{N}(x)dx-\mu\left(\int_{-\infty}^{\infty}\rho_{N}(x)dx-1\right)$
$\displaystyle+\frac{J\text{sgn}(k)N^{2}}{2}PV\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dy\frac{\rho_{N}(x)\rho_{N}(y)}{|x-y|^{k}}.$
(S21)

where $PV$ stands for principal value. As our aim is to find joint
distribution of $(N-1)$ gaps, it will be useful to approximate the principal
value integral in the Eq. (II.4). To do so we assume that there is a particle
at $x$ and break up the $y$ integrals in the three regions,
($-\infty,x-\bar{\Delta}$),($x-\bar{\Delta},x+\bar{\Delta}$) and
($x+\bar{\Delta},\infty$), where $\bar{\Delta}$ is the separation between two
adjacent particles. Noting that the $\rho_{N}(y)$ vanishes in the window
($x-\bar{\Delta},x+\bar{\Delta}$), we rewrite Eq. (II.4) as

$\displaystyle\mathcal{E}[\rho_{N}(x)]\approx\frac{N}{2}\int_{-\infty}^{\infty}x^{2}\rho_{N}(x)dx-\mu\left(\int_{-\infty}^{\infty}\rho_{N}(x)dx-1\right)$
$\displaystyle+\frac{J\text{sgn}(k)N^{2}}{2}\int_{-\infty}^{\infty}dx\Bigg{[}\int_{-\infty}^{x-\bar{\Delta}}+\int_{x+\bar{\Delta}}^{\infty}\Bigg{]}dy\frac{\rho_{N}(x)\rho_{N}(y)}{|x-y|^{k}}.$
(S22)

Figure S2: Schematic description of the discretisation of the integral in the
action $\delta\Sigma$.

The separation $\bar{\Delta}$ at $x$ can be expressed in terms of the density
as

$\bar{\Delta}=\frac{1}{N\rho_{N}(x)}.$ (S23)

Figure S3: Comparison of actual values of the variance of mid-gap calculated
using all three techniques, MC, MH and FT for $\beta=1$.

As the density fluctuates, the gap profile also fluctuates and we are
interested to compute the variance of the gap in the bulk of the system. The
average value of the gap is given by the equilibrium density profile
$\rho_{N}^{\rm(eq)}(x)$ i.e.
$\langle\bar{\Delta}\rangle=\frac{1}{N\rho_{N}^{\rm(eq)}(x)}$. The probability
of a given density profile $\rho_{N}(x)$ Kumar  _et al._ (2020) is

$\mathbb{P}[\rho_{N}]\sim
e^{-\beta\delta\Sigma},~{}\text{with,}~{}\delta\Sigma=\Sigma[\rho_{N}]-\Sigma[\rho_{N}^{\rm(eq)}]$
(S24)

where $\Sigma[\rho_{N}]=\mathcal{E}[\rho_{N}(x)]-\beta^{-1}S[\rho_{N}(x)]$ is
the free energy functional. Writing the fluctuations of density profile about
the equilibrium as $\rho_{N}(x)=\rho^{\rm(eq)}_{N}(x)+\delta\rho(x)$ and using
the relation in Eq. (S23) we obtain

$\bar{\Delta}(x)=\langle\bar{\Delta}\rangle+\delta\Delta=\frac{1}{N\rho_{N}^{\rm(eq)}(x)}-\frac{\delta\rho(x)}{N\left(\rho_{N}^{\rm(eq)}(x)\right)^{2}}.$
(S25)

Using this Eq. (S25) we expand the exponent
$\delta\Sigma=\Sigma[\rho_{N}]-\Sigma[\rho_{N}^{\rm(eq)}]$ in Eq. (S24) in
powers of $\delta\rho(x)$ (up to second order in $\delta\rho$) to get

$\displaystyle\begin{split}\delta\Sigma[\delta\rho,\rho_{N}^{\rm(eq)}]&\approx\frac{J\zeta(k)k(k+1)N^{k+1}}{2}\\\
&\times\int_{-l_{k}}^{l_{k}}dx~{}[\rho_{N}^{\rm(eq)}]^{k-1}~{}\delta\rho(x)^{2}\\\
&+\frac{N}{2\beta}\int_{-l_{k}}^{l_{k}}dx~{}\frac{\delta\rho(x)^{2}}{\rho_{N}^{\rm(eq)}(x)},~{}\quad\text{for}~{}k>1,\end{split}$
(S26)

and

$\displaystyle\begin{split}&\delta\Sigma[\delta\rho,\rho_{N}^{\rm(eq)}]\\\
&\approx\frac{2J\text{sgn}(k)N^{k+1}}{2}\int_{-l_{k}}^{l_{k}}dx~{}[\rho_{N}^{\rm(eq)}]^{k-1}~{}\delta\rho(x)^{2}\\\
&+\frac{J\text{sgn}(k)N^{2}}{2}\int_{-l_{k}}^{l_{k}}dx\left[\int_{-l_{k}}^{x-\langle\bar{\Delta}\rangle}+\int^{l_{k}}_{x+\langle\bar{\Delta}\rangle}\right]dy\\\
&~{}~{}~{}~{}~{}~{}~{}\times~{}\frac{\delta\rho(x)~{}\delta\rho(y)}{|x-y|^{k}},\\\
&~{}~{}+\frac{N}{2\beta}\int_{-l_{k}}^{l_{k}}dx~{}\frac{\delta\rho(x)^{2}}{\rho_{N}^{\rm(eq)}(x)},\quad\text{for}~{}-2<k<1.\end{split}$
(S27)

It is pertinent to note that, in the quadratic approximation of the field
theory, the local term appears due to the fluctations of the coarse grained
gap $\bar{\Delta}$ in the non-local term of the energy
$\mathcal{E}[\rho_{N}(x)]$ for $-2<k<1$. However for $k>1$ the field theory in
Eq. (S20) is local to leading order in $N$. This leads to Eq. (S26) which
further implies $\sigma^{2}_{\Delta_{N/2}}\sim N^{-2}$. We have neglected the
contribution from non-local terms in the field theory because they are
expected to contribute to $\sigma^{2}_{\Delta_{N/2}}$ at most at the same
order ($N^{-2}$) or less, thereby leaving the scaling exponent value $b_{k}$
unchanged. For $N$ (large but finite) particles there are $(N-1)$ number of
gap variables. In order to find the joint distribution of these (discrete) gap
variables from the field theory description, we need to discretise
$\delta\Sigma$ given in Eqs. (S26) and (S27). To do so, we discretise the
integral in the action $\delta\Sigma$ about the equilibrium positions
$\\{y_{i}\\}$ (see Fig. S2). Recall that the microscopic Hessian was computed
about this position configuration in Eq. (7) of the main text. Note that the
equilibrium (in this case the minimum energy) position configuration
$\\{y_{i}\\}$ correspond to equilibrium macroscopic density
$\rho_{N}^{\rm(eq)}(x)$. Hence replacing the integrals as
$\int_{-l_{k}}^{l_{k}}dx\rightarrow\sum_{i}\frac{1}{N\rho^{\rm(eq)}_{N}(y_{i})}$
and performing some simplifications we get the following joint distribution of
the gap variables $\\{\delta\bar{\Delta}_{i}\\}$:

$\displaystyle\mathcal{P}_{\rm
FT}\left(\\{\delta\bar{\Delta}_{i}\\}\right)\sim
e^{-\frac{\beta}{2}\sum_{i,j=1}^{N}M_{ij}\delta\bar{\Delta}_{i}\delta\bar{\Delta}_{j}},\text{where}$
(S28) $\displaystyle
M_{ii}=\begin{cases}J\zeta(k)k(k+1)N^{k+2}[\rho^{\rm(eq)}_{N}(y_{i})]^{k+2}&\text{for
}k>1\\\ 2JN^{k+2}[\rho^{\rm(eq)}_{N}(y_{i})]^{k+2}&\text{for }0<k<1\\\
N^{2}\beta^{-1}[\rho^{\rm(eq)}_{N}(y_{i})]^{2}&\text{for }-2<k<0,\end{cases}$
$\displaystyle M_{ij}=\begin{cases}0&\text{for }k>1\\\
JN^{2}~{}\text{sgn}(k)\frac{\rho^{\rm(eq)}_{N}(y_{i})\rho^{\rm(eq)}_{N}(y_{j})}{|y_{i}-y_{j}|^{k}}&\text{for
}-2<k<1,\end{cases}$ (S29)

in the leading orders in $N$. Using the scaling form
$\rho^{\rm(eq)}_{N}(x)=\frac{1}{l_{k}N^{\alpha_{k}}}F_{k}\left(\frac{x}{l_{k}N^{\alpha_{k}}}\right)$
with $\alpha_{k}=\frac{k}{k+2}$ for $k>1$ and $\frac{1}{k+2}$ for $-2<k<1$, we
find the following system size dependence of the matrix elements (see the next
subsection for details):

$\displaystyle M_{ii}$ $\displaystyle\sim$
$\displaystyle\begin{cases}\mathcal{O}\left(N^{2}\right)&\text{for }k>1\\\
\mathcal{O}\left(N^{k+1}\right)&\text{for }0<k<1\\\
\mathcal{O}\left(N^{2(k+1)/(k+2)}\right)&\text{for }-2<k<0,\end{cases}$ (S30)
$\displaystyle M_{ij}$ $\displaystyle\sim$
$\displaystyle\begin{cases}{0}&\text{for }k>1\\\
\mathcal{O}\left(N^{k+1}\right)&\text{for }-2<k<1,\end{cases}$ (S31)

where we have ignored terms with $|i-j|\sim\mathcal{O}(N)$. From Eq. (S28), it
is easy to compute the variance of middle gap $\Delta_{N/2}$ given by
$\langle\delta\bar{\Delta}_{N/2}^{2}\rangle={[(\beta
M)^{-1}]_{\frac{N}{2}\frac{N}{2}}}$. This can be easily evaluated by
numerically inverting the matrix $M$ and we find the following scaling of the
variance of middle gap

$\displaystyle\sigma^{2}_{\Delta_{N/2}}\sim\begin{cases}\mathcal{O}\left(N^{-2}\right)&\text{for
}k>1\\\ \mathcal{O}\left(N^{-(k+1)}\right)&\text{for }0<k<1\\\
\mathcal{O}\left(N^{-2(k+1)/(k+2)}\right)&\text{for }-1<k<0\\\
\mathcal{O}\left(N^{-(k+1)}\right)&\text{for }-2<k<-1,\end{cases}$ (S32)

which match remarkably well with the large-N scaling obtained from MC
simulation. We notice that the inverse of the dominant term in the matrix
$M_{ij}$ irrespective of whether it is diagonal or off-diagonal gives us the
expected conjectured form for the scaling exponent $b_{k}$.

###  II.5 System size scaling of the matrix elements $M$ and $H$

The variance of the bulk gap and the coarse grained bulk gap are described by
MH and FT respectively and these are in turn governed by the matrices $H$ and
$M$ respectively. Since they quantify the similar physical quantities it is
natural to compare the matrices element wise which is shown in Table (1).

| $M_{ij}$ | $H_{ij}$
---|---|---
Range of $k$ | $i=j$ | $|i-j|\sim\mathcal{O}(1)$ | $i=j$ | $|i-j|\sim\mathcal{O}(1)$
$k>1$ | $N^{2}$ | ${0}$ | $N^{2}$ | $N^{2}$
$0<k<1$ | $N^{1+k}$ | $N^{1+k}$ | $N^{1+k}$ | $N^{1+k}$
$-1<k<0$ | $N^{\frac{2(1+k)}{(k+2)}}$ | $N^{1+k}$ | $N$ | $N^{1+k}$
$-2<k<-1$ | $N^{\frac{2(1+k)}{(k+2)}}$ | $N^{1+k}$ | $N$ | $N^{1+k}$
Table 1: System size dependence of the matrix elements $M_{ij}$ and $H_{ij}$
for different values of $i$ and $j$. As we are interested mainly in the bulk
part of the system we ignore such terms where $(i-j)\sim\mathcal{O}(N)$.

We notice that the elements of matrices $H$ and $M$ have similar system size
scaling for $0<k<1$, while for $k>1$ off-diagonal terms and for $-2<k<0$ the
diagonal terms differ. The variance of the bulk gap are found by the numerical
inversion of these matrices $H$ and $M$. The scaling behaviour of the elements
of the inverse matrices are presented in Table. (2)

| $M^{-1}_{ij}$ | $H^{-1}_{ij}$
---|---|---
Range of $k$ | $i=j$ | $|i-j|\sim\mathcal{O}(1)$ | $i=j$ | $|i-j|\sim\mathcal{O}(1)$
$k>1$ | $N^{-2}$ | ${0}$ | $N^{-2}$ | $N^{-2}$
$0<k<1$ | $N^{-(1+k)}$ | $N^{-(1+k)}$ | $N^{-(1+k)}$ | $N^{-(1+k)}$
$-1<k<0$ | $N^{\frac{-2(1+k)}{(k+2)}}$ | $N^{\frac{-2(1+k)}{(k+2)}}$ | $N^{-(1+k)}$ | $N^{-(1+k)}$
$-2<k<-1$ | $N^{-(1+k)}$ | $N^{-(1+k)}$ | $N^{-(1+k)}$ | $N^{-(1+k)}$

Table 2: System size dependence of the matrix elements $M^{-1}_{ij}$ and
$H^{-1}_{ij}$ for different values of $i$ and $j$. Again we ignore such terms
where $(i-j)\sim\mathcal{O}(N)$.

We observe that the inverse of the dominant term in the matrix $M$ gives the
scaling similar to the variance calculated by numerical inversion of these
matrices. This leads us to the conjecture for the variance given in Eq. (S8).
It is clear from the Table. 2 that except in the regime $-1<k<0$ the
microscopic Hessian theory captures the system size dependence of variance of
bulk gap whereas the quadratic field theory gives correct exponent $b_{k}$
(consistent with MC simulation) for all $k>-2$.

The $N$ dependence of the variance as given in Eq. (S8) is verified
numerically and presented in the Fig. 3 of the main text. However, this does
not contain information on the actual values obtained using different methods
(MC, MH and FT). This comparison is presented in Fig. S3. While the slopes are
consistent with the conjecture, the actual values do not always match which is
not surprising. We find that whenever MH and MC system size scaling agree
[$k\notin(-1,0)$], they even match quantitatively for large-$N$. Remarkably,
in the region $-1<k<0$, the MC and FT results match quantitatively as well.

Figure S4:  Relative fluctuations
$\frac{\sigma_{\Delta_{i}}}{\langle\Delta_{i}\rangle}$ are shown as a function
of $i$ for different system sizes $N$ as well as for different $k$ in the
regime $-1\leq k\leq 0$ for two different system sizes ($N=256$ in pink solid
line and $N=512$ in blue dashed line). For -$1<k<0$, from our numerics we find
that $a_{k}-b_{k}/2$ is slightly different from zero, possibly due to finite
size effects, hence we get a system size dependence in the relative
fluctuation as well. (Inverse temperature $\beta$ is taken to be unity.)

###  II.6 Error estimation in the exponent $b_{k}$

Note that $b_{k}$ is extracted from variance data for last two largest system
sizes. Here we perform error analysis of the exponent value $b_{k}$ for
different $k$ in Fig. 3 of the main text. To do that we first express the
exponent $b_{k}$ in terms of the variance
$\mathcal{V}(N)=\sigma^{2}_{\Delta_{N/2}}$ for a given system size $N$:

$b_{k}=\frac{\ln(\mathcal{V}({N/2}))-\ln(\mathcal{V}({N}))}{\ln(2)}.$ (S33)

The maximum fractional error in the computation of $b_{k}$ can be written in
terms of the error $\delta\mathcal{V}({N/2})$ in the variance
$\mathcal{V}({N})$ as

$\displaystyle(\delta b_{k})_{max}$
$\displaystyle=\left(\frac{\delta\mathcal{V}({N/2})}{\mathcal{V}({N/2})}+\frac{\delta\mathcal{V}({N})}{\mathcal{V}({N})}\right)\times\frac{1}{\ln(2)}.$
(S34)

The error in the variance is given by

$\displaystyle\delta\mathcal{V}({N})\approx\frac{\sqrt{\langle(\Delta_{N/2}-\langle\Delta_{N/2}\rangle)^{4}\rangle-\langle(\Delta_{N/2}-\langle\Delta_{N/2}\rangle)^{2}\rangle^{2}}}{\sqrt{R}},$
(S35)

Figure S5: The distributions of the $i$-th normalized gap
$p^{(i)}_{N}(s)=\langle\delta(s-s_{i})\rangle$ are shown in the regime $-1\leq
k\leq 0$. It is clear that $p^{\rm(i)}_{N}(s)$ is independent of $i$. In these
figures we have taken $N=256$ at inverse temperature $\beta=1$.

where $R$ is the number of independent realisations. We choose $R$ independent
realisations in the simulation by leaving out some Monte-Carlo cycles between
successive data collections so as to ensure negligible correlation. We then
use the above expressions to compute the error bars in Figs. 2 and 3 in the
main text.

##  III Interesting aspects of the gap statistics for $-1\leq k\leq 0$

The system size dependence of relative fluctuations
$\frac{\sigma_{\Delta_{i}}}{\langle\Delta_{i}\rangle}$ in the bulk is
characterized by the exponent $(a_{k}-b_{k}/2)$. From the conjecture Eq. (S8)
the exponent $(a_{k}-b_{k}/2)$ is zero throughout the range $-1\leq k\leq 0$.
We believe that the slight deviations from the predictions for few values of
$k$ are due to finite size effects. Interestingly, in the range $-1\leq k\leq
0$, we find that the relative fluctuations have a weak dependence on $i$, in
the bulk (with almost independent for $k=-1,0^{+}$. This is shown in Fig: S4.
We find that, not only do the relative fluctuations
$\frac{\sigma_{\Delta_{i}}}{\langle\Delta_{i}\rangle}$ have small variations
with $i$, but also the spacing distributions
$p^{\rm(i)}_{N}(s)=\langle\delta(s-s_{i})\rangle$ for individual gaps are also
almost the same (see Fig. S5). However we expect the differences between
$P_{N}(s)$ and $P^{i}_{N}(s)$, for $i$ in bulk, to show up at very small and
very large $s$ where edge statistics (e.g $i=1$) could start dominating over
bulk behaviour.

Figure S6:  Comparison of the level spacing distribution using ‘picket fence
(PF)’ approximation and the actual distribution using Monte-Carlo simulation
(MC). The distribution using PF approximation is fitted with the Eq. (S36)
(Fit) with parameters $A=410.354,B=1.160,C=4.441$. We observe that there is no
convergence of the actual level spacing distribution (MC) with system size. (
Inverse temperature $\beta$ is taken to be unity.)

##  IV Calogero-Moser system

For the CM model (k = 2) it was claimed Bogomolny  _et al._ (2009) that the
spacing distribution follows a form analogous to the WS, namely,

$P(s)=Ae^{-B^{2}/s^{2}-Cs},$ (S36)

where $B=1.46$ is $a$ fitting parameter while $A$ and $C$ are fixed from the
two equations:

$\int_{0}^{\infty}P(s)ds=1\text{~{}and~{}}\int_{0}^{\infty}sP(s)ds=1.$

However, our results for this system differ from this claim. For the CM model,
we found that the distribution, $P_{N}(s)$, of $s$ does not in fact converge
with increasing $N$, which thus indicates that the generalization of the
Wigner surmise does not work here. We believe that the disagreement can be
attributed to the ‘picket fence’ (PF) approximation ($x_{j}\sim j$) used in
the numerical computation of the level spacing distribution in Ref. Bogomolny
_et al._ , 2009. Using this approximation, the Lax matrix $\tilde{L}$ takes
the form

$\tilde{L}_{nm}=p_{n}\delta_{nm}+\frac{i}{2}\frac{1-\delta_{nm}}{(n-m)},$
(S37)

where $p_{n}$ is chosen from a uniform distribution between $-1$ and $1$. From
the eigenvalues of $\tilde{L}$, the level spacing distribution was obtained
and shown to follow the expression Eq. S36. Using this method we can indeed
verify (see Fig. S6) that the level spacing distribution converges to an
$N$-indepependent form that is described by a WS-like distribution. On
comparing the level spacing distribution obtained using the above
approximation (denoted by PF in Fig. S6) with the actual level spacing
distributions (denoted by MC), we observe that they are completely different.
In particular, for the true distribution, we find that there is no convergence
as the system size grows.

## References

  * Bogomolny  _et al._ (2009) E. Bogomolny, O. Giraud, and C. Schmit, Phys. Rev. Lett. 103, 054103 (2009).
  * Riesz (1938) M. Riesz, Acta Sci. Math. Univ. 9, 1 (1938).
  * Agarwal  _et al._ (2019) S. Agarwal, A. Dhar, M. Kulkarni, A. Kundu, S. N. Majumdar, D. Mukamel, and G. Schehr, Phys. Rev. Lett. 123, 100603 (2019).
  * Fletcher (1970) R. Fletcher, The Computer Journal 13, 317 (1970), https://academic.oup.com/comjnl/article-pdf/13/3/317/988678/130317.pdf .
  * Broyden (1970) C. G. Broyden, IMA Journal of Applied Mathematics 6, 76 (1970), https://academic.oup.com/imamat/article-pdf/6/1/76/2233756/6-1-76.pdf .
  * Kumar  _et al._ (2020) A. Kumar, M. Kulkarni, and A. Kundu, Physical Review E 102 (2020), 10.1103/physreve.102.032128.