# Singularities in Hessian element distributions of amorphous media

Vishnu V. Krishnan *  D<EMAIL_ADDRESS>Smarajit Karmakar *  D
<EMAIL_ADDRESS>Kabir Ramola *  D<EMAIL_ADDRESS>Centre for
Interdisciplinary Sciences, Tata Institute of Fundamental Research, Hyderabad
500046, India

###### Abstract

We show that the distribution of elements $H$ in the Hessian matrices
associated with amorphous materials exhibit singularities $P(H)\sim{\lvert
H\rvert}^{\gamma}$ with an exponent $\gamma<0$, as $\lvert H\rvert\to 0$. We
exploit the rotational invariance of the underlying disorder in amorphous
structures to derive these exponents exactly for systems interacting via
radially symmetric potentials. We show that $\gamma$ depends only on the
degree of smoothness $n$ of the potential of interaction between the
constituent particles at the cut-off distance, independent of the details of
interaction in both two and three dimensions. We verify our predictions with
numerical simulations of models of structural glass formers. Finally, we show
that such singularities affect the stability of amorphous solids, through the
distributions of the minimum eigenvalue of the Hessian matrix.

Glass, Hessian

Introduction: Understanding and modelling the properties of amorphous solids
such as glasses has remained a challenge due to their extreme non-equilibrium
nature as well as the underlying disorder in the arrangement of particles [1,
2, *kirkpatrick1988comparison, *kirkpatrick1988mean, *kirkpatrick1989random,
6, 7, 8, 9, 10, 11, 12, 13]. An important aspect in the study of such systems
is their vibrational properties, typically probed through the eigenvalues of
the Hessian matrix [14, 15, 16], with low-frequency modes being particularly
relevant at the low-temperatures where glass physics dominates [17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 16, 31]. Since glass forming
systems settle into disordered configurations, random matrix treatments
provide a natural framework with which to model amorphous systems [32, 33, 34,
35, 36, 37, 38, 39, 40, 41, 42]. The Hessian is thus naturally characterised
by the distribution of its elements [43, 44]. Although the Hessian matrix is
relatively simple to compute in simulations [45, 46], an experimental
determination remains difficult, accessible for example, only through
displacement correlations in colloidal glasses [47]. In this context, it is
important to study the vibrational properties of model systems via
simulations. Simulation models vary in strength and range of interaction [48],
and exhibit a wide range of physical properties [49]. Even though there are
indications of an underlying “amorphous order” in such systems [12, 13], their
behaviour differs fundamentally from that of crystals, with the low lying
eigenvalues of their Hessian matrices displaying marked non-Debye behaviour
[24, 25, 26, 27, 28, 29, 30, 31], related to the growth of a structural length
scale [50] as well as static correlations [51]. The distribution of elements
affects the eigenvalues of the Hessian which is directly related to the
mechanical properties of amorphous solids [52, 53, 54].

(a) Crystal$-$250$-$200$-$150$-$100$-$50050$\displaystyle H$(b)
Amorphous$\displaystyle P\left(H\right)$

Figure 1: Comparison of energy minimised configurations in (a) a mono-disperse
crystal and (b) an amorphous glass consisting of two types of particles: A
(small, green) and B (large, blue). Also plotted are the corresponding
distributions $P(H)$ of the diagonal elements of their Hessians. In (b), the
distribution is shown only for A-B interactions.

In this Letter, we analyse the distributions of Hessian elements in structural
glass formers analytically, as well as through numerical simulations. In
simulations, interaction potentials are typically cut-off at a finite distance
for computational expediency, and in order to avoid unphysical changes, they
are smoothed to relevant degrees at this cut-off [55]. However, the effect of
this smoothness on vibrational properties of such amorphous systems has never
been investigated. In crystals, the ordered inter-particle distances imply
delta distributed Hessian elements. For example, in a triangular lattice of
particles with nearest neighbour interactions as shown in Fig. 1 (a), the
diagonal Hessian elements take on one of two non-zero values. However, in the
case of a disordered arrangement derived from simulations of glass formers as
shown in Fig. 1 (b), this distribution is continuous, peaking at zero with a
marked singularity. We exploit the rotational invariance of the underlying
amorphous disorder to derive analytic expressions for these distributions, for
systems with radially symmetric interactions, which we then verify using
direct numerical simulations. We show that these distributions indeed exhibit
smoothness-dependent singularities

$P\left(H\right)\sim{\left\lvert
H\right\rvert}^{\gamma},~{}~{}~{}\text{for}~{}~{}~{}|H|\rightarrow 0$ (1)

with a power $\gamma<0$. We derive exact results for $\gamma$, for any degree
of smoothness, in both two dimensions (2D) as well as three dimensions (3D).
Our results are summarized in Table 1, highlighting the non-trivial dependence
of these distributions on the nature of the interaction at the cut-off.
Finally, we show that these singularities have crucial implications for the
low-energy vibrational modes of amorphous systems, through a numerical
sampling of the minimum eigenvalue of systems with varying smoothness in their
interactions.

Distribution of Hessian Elements: We begin with a disordered arrangement of
particles at positions $\\{\mathbf{r}^{i}\\}$, where the particle index
$i\in\\{1,\ldots,N\\}$. The Hessian is then a $dN\times dN$ matrix where $d$
represents the dimension of the system, with elements
$\mathcal{H}^{ij}_{\alpha\beta}$ that describe the stiffness between particles
$i$ and $j$, along the coordinates $\alpha,\beta\in\\{x,y,z\\}$. In terms of
the inter-particle distance vector
$\mathbf{r}^{ij}=\mathbf{r}^{i}-\mathbf{r}^{j}$, the Hessian elements may be
expressed as

$\mathcal{H}_{\alpha\beta}^{ij}(\mathbf{r}^{ij})=\frac{\partial^{2}U\left[\\{\mathbf{r}^{i}\\}\right]}{\partial
r^{ij}_{\alpha}\partial r^{ij}_{\beta}}.$ (2)

Here, $U\left[\\{\mathbf{r}^{i}\\}\right]$ is the total potential energy of
the system, which is a function of the positions of all particles. For
pairwise additive interactions,
$U\left[\\{\mathbf{r}^{ij}\\}\right]=\sum_{ij}\psi^{ij}$, where $\psi^{ij}$ is
the interaction potential between particles $i$ and $j$. We consider central
potentials $\psi(r)$, where $r\equiv r^{ij}=\lvert\mathbf{r}^{ij}\rvert$ is
the distance between particles $i\text{ and }j$. In addition, these potentials
are smoothed to $n$ derivatives at a cut-off distance $r_{c}$, i.e.,
$\left.\frac{d^{m}\psi}{dr^{m}}\right\rvert_{r_{c}}=0$ for all $0\leq m\leq
n$. Fig. 2 shows a typical interaction potential in 2D used in our numerical
simulations, and its derivatives, all of which tend to zero at cut-off ($n=2$
in this case). The Hessian elements for central potentials are given by [56]

$\mathcal{H}_{\alpha\beta}^{ij}(\mathbf{r}^{ij})=-\left(\frac{\psi_{rr}^{ij}}{{(r^{ij})}^{2}}-\frac{\psi_{r}^{ij}}{{(r^{ij})}^{3}}\right)r^{ij}_{\alpha}r^{ij}_{\beta}-\delta_{\alpha\beta}\frac{\psi_{r}^{ij}}{r^{ij}},$
(3)

where the subscripts of $r$ indicate partial derivatives with respect to the
inter-particle distances:
$\psi_{r}\equiv\psi_{r}^{ij}=\partial_{r^{ij}}\psi^{ij},\;\psi_{rr}\equiv\psi_{rr}^{ij}=\partial_{r^{ij}}\partial_{r^{ij}}\psi^{ij}$.
Since the potential and its derivatives vanish at the cut-off distance, the
small Hessian elements arise primarily due to pair-distances near this cut-
off. The values of Hessian elements in Eq. (3) depend on the length and angle
of the inter-particle distances in an arbitrarily chosen (fixed) Cartesian
coordinate system. We therefore define a generalised angular coordinate
$\Omega$ with respect to these fixed axes, which is a function of the angle
$\phi$ in 2D, and both the polar and azimuthal angles $\theta,\phi$ in 3D. The
distribution of Hessian elements $P(H)$ is then given by

$P\left(H\right)=\int
dr\;d\Omega\;P\left(r,\Omega\right)\delta{\left(H-\mathcal{H}_{\alpha\beta}^{ij}(\mathbf{r}^{ij})\right)},$
(4)

where $P(r,\Omega)$ represents the joint probability distributions of inter-
particle distances and orientations, as shown in the inset of Fig. 3. For a
given $H$, the delta function constraint in Eq. (4) selects the corresponding
contour in $(r,\Omega)$ through Eq. (3), as shown in the inset of Fig. 2.

0.81.01.21.41.6$\displaystyle
r$0100200300400$\displaystyle\psi\left(r\right)$CUT-
OFF$\displaystyle(r_{c})$$\displaystyle\psi$$\displaystyle
10^{-2}\psi_{r}$$\displaystyle
10^{-2}\psi_{rr}$$-$1.0$-$0.50.00.51.0$\displaystyle\cos\left(\phi\right)$1.01.21.41.6$\displaystyle
r$$\displaystyle H=0$$\displaystyle H>0$$\displaystyle H<0$$\displaystyle
r_{c}$$\displaystyle H_{\alpha\alpha}$$-$25$-$20$-$15$-$10$-$50510

Figure 2: The R10 potential $\psi(r)\sim r^{-10}+c+br^{2}+ar^{4}$ for A-B
interactions, smoothed to two derivatives $(n=2)$ at the cut-off
$[\psi(r_{c})=\psi^{\prime}(r_{c})=\psi^{\prime\prime}(r_{c})=0]$. The
derivatives have been scaled down by a factor of 100 for clarity. The inset
shows a plot in the inter-particle coordinates $(r,\cos\phi)$, displaying
contours which contribute to the distribution of diagonal Hessian elements
($\alpha=\beta$) at a fixed $H$.

We next exploit the rotational invariance of the underlying amorphous disorder
to relate $P(H)$ to the distribution of inter-particle distances. This
microscopic rotational symmetry leads to isotropic angular distributions at
all distances. Hence, we may assume the distributions of inter-particle
distances and angles to be uncorrelated $P(r,\Omega)=P(r)P(\Omega)$. Numerical
sampling of this joint distribution $P(r,\Omega)$ in structural glass formers
(inset of Fig. 3), demonstrates that all orientations are sampled uniformly,
independent of the radial distance. We therefore assume uniform distributions
in the angular variables
$P_{2\text{D}}(\phi)=P_{3\text{D}}(\phi)=\frac{1}{2\pi}$, and
$P_{3\text{D}}(\cos\theta)=\frac{1}{2}$. The distribution of $r$ is the radial
distribution function $g_{2}(r)$, normalised over the range of interaction
$r\in[0,r_{c}]$, i.e., $P(r)=g_{2}(r)\Theta(r_{c}-r)\mathcal{N}$, where
$\Theta$ is the Heaviside function, and
${\mathcal{N}}^{-1}=\int_{0}^{r_{c}}g_{2}(r)$ is a constant of normalisation.
For a given $H$, the integration over the angular variables in Eq. (4) yields
a Jacobian factor, leading to the Hessian element distribution

$P\left(H\right)=\int_{0}^{r_{c}}dr\;P(r)\;\frac{P\left(\Omega\right)}{\left\lvert\frac{\partial
H}{\partial\Omega}\right\rvert}=\int_{0}^{r_{c}}dr\;P(r)\;\mathcal{P}(H,r).$
(5)

Above, the angular variable $\Omega$ in the final expression has been
expressed in terms of $H$ and $r$ using Eq. (3). Therefore, given an empirical
$g_{2}(r)$, $P(H)$ can be determined exactly. The partial integrand
$\mathcal{P}(H,r)={\left[\left\lvert\frac{\partial
H}{\partial\Omega}\right\rvert\right]}^{-1}P\left(\Omega\right)$ depends on
the dimension, as well as the spatial indices $\alpha$ and $\beta$. However,
rotational symmetry dictates that there are only two classes of Hessian
element distributions: the diagonal ($H_{\alpha\alpha}:H_{xx}\equiv
H_{yy}\equiv H_{zz}$) and off-diagonal ($H_{\alpha\beta}:H_{xy}\equiv
H_{yz}\equiv H_{zx}$).

Exact Partial Integrands: We first consider the distribution of Hessian
elements in 2D systems, choosing the diagonal element $H_{xx}$ and the off-
diagonal element $H_{xy}$. From Eq. (3),
$H^{2\text{D}}_{xx}=-\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\cos^{2}\phi-\frac{\psi_{r}}{r}$,
and
$H^{2\text{D}}_{xy}=-\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\cos\phi\sin\phi$.
It is convenient to use the angular variable $\Omega\equiv\cos\phi$ with the
distribution
$P_{2\text{D}}\left(\Omega\right)=1/\left(\pi\sqrt{1-\Omega^{2}}\right)$.
Substitution of these into the expression in Eq. (5) yields the same partial
integrand for both cases
$\mathcal{P}^{2\text{D}}_{\alpha\beta}={\left[\pi\left\lvert\frac{\partial
H}{\partial\cos\phi}\right\rvert\sqrt{1-\cos^{2}\phi}\right]}^{-1}$. Finally,
inverting the above expressions for $H^{2\text{D}}_{xx}$ and
$H^{2\text{D}}_{xy}$ to express $\cos\phi$ in terms of $H$ and $r$, we arrive
at

$\displaystyle\mathcal{P}^{2\text{D}}_{\alpha\beta}\left(H,r\right)=\begin{cases}{\left\lvert
4\pi^{2}\left(H+\frac{\psi_{r}}{r}\right)\left(H+\psi_{rr}\right)\right\rvert}^{-1/2}&\alpha=\beta,\\\
{\left\lvert\pi^{2}\left\\{{\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}^{2}-4H^{2}\right\\}\right\rvert}^{-1/2}&\alpha\neq\beta.\end{cases}$
(6)

For 3D systems, we choose $H_{zz}$ and $H_{xz}$ to represent the diagonal and
off-diagonal elements, respectively. From Eq. (3) we have
$H^{3\text{D}}_{zz}=-\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\cos^{2}\phi-\frac{\psi_{r}}{r}$,
and
$H^{3\text{D}}_{xz}=-\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\cos\theta\sin\theta\cos\phi$.
In the case of diagonal elements, it is convenient to choose
$\Omega\equiv\cos\theta$, with the distribution
$P_{3\text{D}}\left(\cos\theta\right)=\frac{1}{2}$ (see inset of Fig. 3).
Under these conditions, the partial integrand defined in Eq. (5) is
$\mathcal{P}^{3\text{D}}_{zz}={\left[2\left\lvert\frac{\partial
H}{\partial\cos\theta}\right\rvert\right]}^{-1}$. In the case of off-diagonal
elements, it is convenient to choose $\Omega\equiv\cos\phi$, and given that
$P\left(\phi\right)=\frac{1}{2}$, we arrive at the exact integral form
$\mathcal{P}^{3\text{D}}_{xz}=\int\pi
d(\cos\theta){\left[\left\lvert\frac{\partial
H}{\partial\cos\phi}\right\rvert\sqrt{1-\cos^{2}\phi}\right]}^{-1}$. Finally,
inverting the above expressions for $H^{3\text{D}}_{zz}\text{ and
}H^{3\text{D}}_{xz}$ to express $\cos\theta\text{ and }\cos\phi$ in terms of
$H$ and $r$, we arrive at the simplified expressions

$\displaystyle\mathcal{P}^{3\text{D}}_{\alpha\beta}\left(H,r\right)=\begin{cases}{\left\lvert
4\left(H+\frac{\psi_{r}}{r}\right)\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\right\rvert}^{-1/2}&\alpha=\beta,\\\
\frac{\kappa}{H}\int^{1}_{-1}dx{\left[x^{2}(1-x^{2})-\kappa^{2}\right]}^{-1/2}&\alpha\neq\beta.\\\
\end{cases}$ (7)

with $\kappa=H{\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}^{-1}$. The above
integral form for the off-diagonal elements has the asymptotic behaviour
(refer to Supplemental Material [57]/.1)

$\mathcal{P}^{3\text{D}}_{\alpha\neq\beta}\left(H,r\right)\overset{\kappa\to
0}{\sim}\frac{2\kappa}{H}\log{\left(\frac{4}{\kappa}\right)}.$ (8)

012345$\displaystyle r$01234567$\displaystyle
g_{2}\left(r\right)$$\displaystyle r_{c}$CUT-OFF\-
$\displaystyle\pi$0$\displaystyle\pi$$\displaystyle\phi$1.251.501.752.002.252.50$\displaystyle
r$$\displaystyle P(r,\phi)$$-$101$\displaystyle\cos(\theta)$$\displaystyle
P(r,\cos\theta)$0.00.20.40.60.81.0

Figure 3: Numerically sampled radial distribution function of A-B particle
pairs in the R10 system in three dimensions. The vertical line indicates the
cut-off distance of the interaction. Inset: Numerically sampled joint
distributions of the distances ($r$) and angles ($\phi$ and $\cos\theta$),
indicating isotropic angular distributions.

The distributions of Hessian elements can now be found using these partial
integrands [Eqs. (6) and (7)] in Eq. (5) along with the numerically obtained
radial distribution functions, as shown in Fig. 3. We display the match
between our theoretical predictions and distributions obtained from numerical
simulations in Fig. 4. In the Supplemental Material [57]/.2, we additionally
show a precise match between $P(H)$ obtained analytically and numerically, for
an exponentially decaying $g_{2}(r)$.

Asymptotic Forms and Singular Distributions: We focus next on the behaviour of
$P(H)$ in the limit $H\to 0$. It is clear from the forms of the partial
integrands in Eqs. (6) and (7), that they diverge as $r\to r_{c}$ and $H\to
0$, since $\psi_{r}\to 0~{}\text{ and }~{}\psi_{rr}\to 0$. Consequently, the
corresponding integrals determining $P(H)$ in Eq. (5) exhibit singularities
$P(H)\sim H^{\gamma}$ as $H\to 0$, as displayed in Figs. 1 (b) and 4. We show
below that the strength $\gamma$ of this singularity depends on the rate of
decay (i.e., smoothness $n$) of the interaction at cut-off.

To determine the behaviour of the distribution near the singularity, we focus
on the integral in Eq. (5), as we approach the cut-off $r\to r_{c}$. It is
assumed also, that the $g_{2}(r)$ does not contain singularities at $r_{c}$ as
justified by the plot in Fig. 3. The interaction potential, smooth to $n$
derivatives at the cut-off distance and its derivatives in the limit of $r\to
r_{c}$, may therefore be approximated as

$\displaystyle\psi\left(r\right)$
$\displaystyle={\left(r_{c}-r\right)}^{n+1}f\left(r\right),$
$\displaystyle\psi_{r}/r$ $\displaystyle\approx
C_{1}{\left(r_{c}-r\right)}^{n},$ $\displaystyle\psi_{rr}$
$\displaystyle\approx C_{2}{\left(r_{c}-r\right)}^{n-1},$ (9)

where $f(r)$ is a regular function, along with $C_{1}=-(n+1)f(r)/r$, and
$C_{2}=n(n+1)f(r)$ which vary significantly slower compared to the power-law
term as $r\to r_{c}$. Above, in $\psi_{r}/r$, we have ignored
${(r_{c}-r)}^{n+1}f^{\prime}(r)$ in comparison to ${(r_{c}-r)}^{n}f(r)$. Under
this approximation, the singular points of $\mathcal{P}(H,r)$ in the
complex-$r$ plane are determined by the two expressions

${\left(r_{c}-r\right)}^{n}=-\frac{H}{C_{1}},\quad{\left(r_{c}-r\right)}^{n-1}=-\frac{H}{C_{2}}.$
(10)

The full structure of the poles in Eq. (10) is detailed in the Supplemental
Material [57]/.3. However, as the integral in Eq. (5) is performed over the
real interval $[0,r_{c}]$, the new upper limit is determined by the largest
positive real root of these expressions. This occurs at a value
$r^{*}=r_{c}-s$, where $s\equiv s(H)$ represents the shift in the upper limit.
The singular behaviour in $\mathcal{P}(H,r)$ is thus determined by both the
sign of $H$, and the signs of the first two derivatives of the potential
($C_{1}$ and $C_{2}$), as the cut-off distance $r\to r_{c}$ is approached.

Table 1: Asymptotic behaviour of Hessian element distributions in the limit $H\to 0$. The diagonal element distribution, depends on the relative signs of $H$ and the interaction potential near the cut-off $\psi_{\delta}\equiv\psi(r_{c}-\delta)$. The case $(*)$ corresponds to $H\times\psi_{\delta}>0$, while $(\dagger)$ corresponds to $H\times\psi_{\delta}<0$. The results are identical for both two and three dimensions. Element $(H)$ | Smoothness $(n)$ | $\lim\limits_{H\to 0}P(H)$
---|---|---
Diagonal | ${[2,\infty)}^{*}$ | ${\left\lvert H\right\rvert}^{-1+\frac{3}{2n}}$
$(\alpha=\beta)$ | ${\\{2\\}}^{\dagger}$ | ${\left\lvert H\right\rvert}^{-1+\frac{3}{2n}}$
| ${\\{3\\}}^{\dagger}$ |  ${\left\lvert H\right\rvert}^{-\frac{1}{2}}$ $\log{\left({\left\lvert H\right\rvert}^{-1}\right)}$
| ${(3,\infty)}^{\dagger}$ | ${\left\lvert H\right\rvert}^{-1+\frac{1}{n-1}}$
Off-diagonal | $\\{2\\}$ | $\log{\left({\left\lvert H\right\rvert}^{-1}\right)}$
$(\alpha\neq\beta)$ | $(2,\infty)$ | ${\left\lvert H\right\rvert}^{-1+\frac{1}{n-1}}$

In order to analyse the asymptotic forms of $P(H)$ in Eq. (5) as $H\to 0$, we
define a small variable $\epsilon$ as a distance to $r^{*}$ at a given value
of $H\text{, }\epsilon=(r_{c}-s)-r$. The derivatives of the potential
described in Eq. (9) can then be written as

$\psi_{r}/r\approx C_{1}{\left(\epsilon+s\right)}^{n},\quad\psi_{rr}\approx
C_{2}{\left(\epsilon+s\right)}^{n-1}.$ (11)

We can now extract the asymptotic behaviour of $P(H)$ using this approximation
in the partial integrands in Eqs. (6) and (7). The behaviour of these
integrands depends on the relative signs of $H$ and the interaction potential
near cut-off $\psi_{\delta}\equiv\psi(r_{c}-\delta)$ with $\delta/r_{c}\ll 1$.
Below, we present the analysis for the diagonal elements in 2D, with the
details of all cases presented in the Supplemental Material [57]/.4. Using Eq.
(11) in Eq. (6), for $H\to 0,\;\epsilon\to 0$, we have

$\mathcal{P}^{2\text{D}}_{\alpha\alpha}\left(H,\epsilon\right)\sim{\left\\{\left[H+C_{1}{(\epsilon+s)}^{n}\right]\left[H+C_{2}{(\epsilon+s)}^{n-1}\right]\right\\}}^{-1/2}.$
(12)

We can extract the limiting behaviour of $P(H)$ by identifying the dominant
contribution from the above expression to the integral in Eq. (5). There exist
up to three regimes of $\epsilon$ in terms of the behaviour of the partial
integrand in Eq. (12): (i) $\left[0,H^{1/n-1}\right)$, (ii)
$\left[H^{1/n-1},H^{1/n}\right)$, and (iii) $\left[H^{1/n},\infty\right)$,
depending on the shift $s(H)$ (refer to Supplemental Material [57]/.3). For
the case $H\times\psi_{\delta}>0$, contribution of the integral over the last
interval dominates, with
$P\left(H\right)\sim\int_{H^{1/n}}^{\infty}\epsilon^{-n+1/2}=\left.{\epsilon^{-n+3/2}}\right|_{H^{1/n}}^{\infty}$.
Therefore, we arrive at the asymptotic form
$P^{2\text{D}}_{\alpha\alpha}\left(H\right)\sim H^{-1+\frac{3}{2n}}$. Our
results for all cases are summarized in Table 1. Remarkably, although the
expressions in Eqs. (6) and (7) have very different forms, they yield exactly
the same results for the singularities in both 2D and 3D.

$\displaystyle{-10^{1}}$$\displaystyle{-10^{-2}}$$\displaystyle{-10^{-5}}$$\displaystyle{-10^{-8}}$$\displaystyle{0}$$\displaystyle{10^{-8}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{1}}$$\displaystyle
H$$\displaystyle{10^{-7}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-3}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{1}}$$\displaystyle{10^{3}}$$\displaystyle
P\left(H\right)$$\displaystyle n=$ 2$\displaystyle n=$
3$\displaystyle{\left|H\right|}^{-1/4}$$\displaystyle{\left|H\right|}^{-1/2}$$\displaystyle{10^{-3}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{0}}$$\displaystyle\lambda_{\min}$$\displaystyle{10^{-3}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{0}}$$\displaystyle
P\left(\lambda_{\min}\right)$

Figure 4: Distribution of diagonal Hessian elements in the R10 system with
different smoothness at cut-off $n=2$ and $3$ in two dimensions. The
configurations sampled are at energy minima. The plot is in symmetric-log-log
scale, showing the power law divergence in both positive and negative $H$. The
dashed lines represent the analytic predictions from Eq. (5), using
numerically sampled radial distribution functions. Exponent of the singularity
changes with the degree of smoothness $n$. (Inset) Distribution of the minimum
eigenvalue $\lambda_{\min}$ of the Hessian (for system size $N=256$),
displaying significant changes with the degree of smoothness $n$.

Numerical Simulations: In order to verify our predictions, we have performed
extensive numerical simulations of structural glass formers, in 2D as well as
3D. We simulate a binary mixture of purely repulsive particles of type A and B
(refer to Supplemental Material [57]/.5). The A-B interactions are illustrated
in Fig. 2. High temperature molecular dynamics simulations were utilised to
generate independent, uncorrelated configurations of particles, and sample
their inherent structures by locating the nearest local minimum via the
conjugate-gradient minimization. We then evaluate the Hessian elements between
particles within interacting range of each other. In Fig. 4 we plot the
numerically sampled distributions of the diagonal Hessian elements for A-B
interactions in 2D systems for $n=2$ and $n=3$, along with our theoretical
predictions (dashed lines) using Eq. (5) in Eq. (6), displaying near-perfect
agreement. These distributions diverge as $H\to 0^{+}$ with the exponents
$-\frac{1}{4}$ for $n=2$ and $-\frac{1}{2}$ for $n=3$, as predicted in Table
1.

Finally, we turn our attention to the vibrational properties of the system,
that are probed through the eigenvalue spectrum of the Hessian matrix. We test
the sensitivity of low-lying eigenvalues to the smoothness of the potential,
and consequently, the power of the singularity in $P(H)$, by analysing the
distribution of the minimum eigenvalue, $\lambda_{\min}$. Our results for the
short-ranged R10 glass model, displayed in the inset of Fig. 4, show a
significant divergence at low values, between distributions for two values of
smoothness.

Discussion: We have presented analytic results for the distribution of Hessian
elements in disordered amorphous media in 2D and 3D, and verified them with
extensive numerical simulations. Our treatment is quite general, relying only
on the isotropy of the underlying amorphous medium, and can be extended to
other systems displaying such disorder. Additionally, we have shown that the
Hessian matrices of amorphous materials display a preponderance of small
elements, characterized by a singularity whose strength depends on the
smoothness of the interaction potential at the cut-off distance. Remarkably,
the results for the singularities are exactly the same in both 2D and 3D, a
fact that warrants deeper investigation. We have also shown numerically that
such singularities have crucial implications for the low-lying eigenvalues of
the Hessian matrix that govern the stability or fragility of amorphous solids,
highlighting their sensitivity to the small, non-zero elements in the Hessian
matrix. Our results are particularly relevant for numerical studies of
glasses, where the degree of smoothness in interaction potentials have been
shown to affect vibrational properties [58]. The limit where such interaction
potentials display sharp cut-offs are typically used in the study of jamming
transitions [59, 60, 61, 62], and it would be interesting to study the effect
of smoothness in the potential on the properties of such systems. Finally, it
would also be interesting to extend our analytic results to construct bounds
on the vibrational density of states of amorphous systems [63].

Acknowledgements: We thank Edan Lerner, Prathyush Manchala, Mustansir Barma,
Daan Frenkel, Srikanth Sastry, Chandan Dasgupta, Satya Majumdar, Kedar Damle
and Deepak Dhar for useful discussions. V. V. K. thanks CSIR for the Shyama
Prasad Mukherjee Fellowship (SPM-07/1142(0228)/2015-EMR-1). S. K. would like
to acknowledge the support from Swarna Jayanti Fellowship Grants No.
DST/SJF/PSA-01/2018-19 and No. SB/SFJ/2019-20/05. This project was funded by
intramural funds at TIFR Hyderabad from the Department of Atomic Energy (DAE).

## References

  * Angell [1988] C. A. Angell, J. Phys. Chem. Solids 49, 863 (1988).
  * Kirkpatrick and Thirumalai [1987] T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. B 36, 5388 (1987).
  * kir [1988a] Phys. Rev. A 37, 4439 (1988a).
  * kir [1988b] Phys. Rev. B 37, 5342 (1988b).
  * kir [1989] J. Phys. A: Mathematical and General 22, L149 (1989).
  * Debenedetti [1996] P. G. Debenedetti,  _Metastable Liquids: Concepts and Principles_ (Princeton University Press, Princeton, NJ, 1996).
  * Das [2004] S. P. Das, Rev. Mod. Phys. 76, 785 (2004).
  * Götze [2008] W. Götze, in  _International Series of Monographs on Physics_, Vol. 143 (Oxford University Press, Oxford, 2008).
  * Shintani and Tanaka [2008] H. Shintani and H. Tanaka, Nature Mater. 7, 870 (2008).
  * Berthier and Biroli [2011] L. Berthier and G. Biroli, Rev. Mod. Phys. 83, 587 (2011).
  * Wolynes and Lubchenko [2012] P. G. Wolynes and V. Lubchenko,  _Structural Glasses and Supercooled Liquids: Theory, Experiment, and Applications_ (Wiley, New York, 2012).
  * Karmakar  _et al._ [2014] S. Karmakar, C. Dasgupta, and S. Sastry, Annu. Rev. Condens. Matter Phys. 5, 255 (2014).
  * Karmakar  _et al._ [2015] S. Karmakar, C. Dasgupta, and S. Sastry, Rep. Prog. Phys. 79, 016601 (2015).
  * Sastry  _et al._ [1997] S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Phys. Rev. E 56, 5533 (1997).
  * Angell  _et al._ [2003] C. A. Angell, Y. Yue, L.-M. Wang, J. R. D. Copley, S. Borick, and S. Mossa, J. Phys.: Condensed Matter 15, S1051 (2003).
  * Zaccone [2020] A. Zaccone, J. Phys.: Condensed Matter 32, 203001 (2020).
  * Angell and Rao [1972] C. A. Angell and K. J. Rao, J. Chem. Phys. 57, 470 (1972).
  * Goldstein [1976] M. Goldstein, J. Chem. Phys. 64, 4767 (1976).
  * Buchenau  _et al._ [1984] U. Buchenau, N. Nücker, and A. J. Dianoux, Phys. Rev. Lett. 53, 2316 (1984).
  * Rufflé  _et al._ [2003] B. Rufflé, M. Foret, E. Courtens, R. Vacher, and G. Monaco, Phys. Rev. Lett. 90, 095502 (2003).
  * Gurevich  _et al._ [2003] V. L. Gurevich, D. A. Parshin, and H. R. Schober, Phys. Rev. B 67, 094203 (2003).
  * Rufflé  _et al._ [2006] B. Rufflé, G. Guimbretière, E. Courtens, R. Vacher, and G. Monaco, Phys. Rev. Lett. 96, 045502 (2006).
  * Moriel  _et al._ [2019] A. Moriel, G. Kapteijns, C. Rainone, J. Zylberg, E. Lerner, and E. Bouchbinder, J. Chem. Phys. 151, 104503 (2019).
  * Monaco and Mossa [2009] G. Monaco and S. Mossa, Proc. Natl. Acad. Sci. U.S.A. 106, 16907 (2009).
  * Monaco and Giordano [2009] G. Monaco and V. M. Giordano, Proc. Natl. Acad. Sci. U.S.A. 106, 3659 (2009).
  * DeGiuli  _et al._ [2014] E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner, and M. Wyart, Soft Matter 10, 5628 (2014).
  * Charbonneau  _et al._ [2016] P. Charbonneau, E. I. Corwin, G. Parisi, A. Poncet, and F. Zamponi, Phys. Rev. Lett. 117, 045503 (2016).
  * Lerner  _et al._ [2016] E. Lerner, G. Düring, and E. Bouchbinder, Phys. Rev. Lett. 117, 035501 (2016).
  * Kapteijns  _et al._ [2018] G. Kapteijns, E. Bouchbinder, and E. Lerner, Phys. Rev. Lett. 121, 055501 (2018).
  * Wang  _et al._ [2019] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, Nat. Commun. 10, 1 (2019).
  * Shimada  _et al._ [2020] M. Shimada, H. Mizuno, L. Berthier, and A. Ikeda, Phys. Rev. E 101, 052906 (2020).
  * Stanifer  _et al._ [2018] E. Stanifer, P. K. Morse, A. A. Middleton, and M. L. Manning, Phys. Rev. E 98, 042908 (2018).
  * Beltukov and Parshin [2011] Y. M. Beltukov and D. A. Parshin, Phys. Solid State 53, 151 (2011).
  * Manning and Liu [2015] M. L. Manning and A. J. Liu, Europhys. Lett. 109, 36002 (2015).
  * Beltukov  _et al._ [2013] Y. M. Beltukov, V. I. Kozub, and D. A. Parshin, Phys. Rev. B 87, 134203 (2013).
  * Baggioli  _et al._ [2019] M. Baggioli, R. Milkus, and A. Zaccone, Phys. Rev. E 100, 062131 (2019).
  * Conyuh  _et al._ [2017] D. A. Conyuh, Y. M. Beltukov, and D. A. Parshin, J. Phys.: Conf. Ser. 929, 012031 (2017).
  * Brito  _et al._ [2010] C. Brito, O. Dauchot, G. Biroli, and J.-P. Bouchaud, Soft Matter 6, 3013 (2010).
  * Altieri [2019] A. Altieri, An exactly solvable model: The perceptron, in  _Jamming and Glass Transitions: In Mean-Field Theories and Beyond_ (Springer, New York, 2019) pp. 65–113.
  * Maradudin and Weiss [1958] A. Maradudin and G. H. Weiss, J. Soc. Indust. Appl. Math. 6, 302 (1958).
  * Gurarie and Chalker [2003] V. Gurarie and J. T. Chalker, Phys. Rev. B 68, 134207 (2003).
  * Benetti  _et al._ [2018] F. P. C. Benetti, G. Parisi, F. Pietracaprina, and G. Sicuro, Phys. Rev. E 97, 062157 (2018).
  * Cavagna  _et al._ [1999] A. Cavagna, I. Giardina, and G. Parisi, Phys. Rev. Lett. 83, 108 (1999).
  * Huang and Wu [2009] B. J. Huang and T.-M. Wu, Phys. Rev. E 79, 041105 (2009).
  * Kittel [2004] C. Kittel,  _Introduction to Solid State Physics_ , 8th ed. (Wiley, New York, 2004).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Cengage Learning, Boston, 1976).
  * Ghosh  _et al._ [2010] A. Ghosh, V. K. Chikkadi, P. Schall, J. Kurchan, and D. Bonn, Phys. Rev. Lett. 104, 248305 (2010).
  * Wang  _et al._ [2020] X. Wang, S. Ramírez-Hinestrosa, J. Dobnikar, and D. Frenkel, Phys. Chem. Chem. Phys. 22, 10624 (2020).
  * Dauchot  _et al._ [2011] O. Dauchot, S. Karmakar, I. Procaccia, and J. Zylberg, Phys. Rev. E 84, 046105 (2011).
  * Karmakar  _et al._ [2012] S. Karmakar, E. Lerner, and I. Procaccia, Physica A 391, 1001 (2012).
  * Rainone  _et al._ [2020] C. Rainone, E. Bouchbinder, and E. Lerner, Proc. Natl. Acad. Sci. U.S.A. 117, 5228 (2020).
  * Maloney and Lemaître [2004] C. Maloney and A. Lemaître, Phys. Rev. Lett. 93, 195501 (2004).
  * Maloney and Lemaître [2006] C. E. Maloney and A. Lemaître, Phys. Rev. E 74, 016118 (2006).
  * Dasgupta  _et al._ [2012] R. Dasgupta, S. Karmakar, and I. Procaccia, Phys. Rev. Lett. 108, 075701 (2012).
  * Frenkel and Smit [2002] D. Frenkel and B. Smit,  _Understanding Molecular Simulation: From Algorithms to Applications_, 2nd ed. (Elsevier, New York, 2002).
  * Karmakar  _et al._ [2010] S. Karmakar, E. Lerner, and I. Procaccia, Phys. Rev. E 82, 026105 (2010), Appendix C, arXiv:1004.2198 [cond-mat.stat-mech] .
  * [57] See Supplemental Material for details, which includes Refs. [60-67].
  * Shimada  _et al._ [2018] M. Shimada, H. Mizuno, and A. Ikeda, Phys. Rev. E 97, 022609 (2018).
  * O’Hern  _et al._ [2003] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
  * Silbert  _et al._ [2006] L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 73, 041304 (2006).
  * Ramola and Chakraborty [2016] K. Ramola and B. Chakraborty, J. Stat. Mech.: Theory Exp. 2016 (11), 114002.
  * Ramola and Chakraborty [2017] K. Ramola and B. Chakraborty, Phys. Rev. Lett. 118, 138001 (2017).
  * Horn and Johnson [2012] R. A. Horn and C. R. Johnson,  _Matrix Analysis_ (Cambridge University Press, Cambridge, UK, 2012) Chap. 6.
  * Plimpton [1995] S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * lam [2003] Large-scale Atomic/Molecular Massively Parallel Simulator, Sandia National Laboratories (2003), Versions: 12Dec2018, 5Jun2019.
  * Anderson  _et al._ [1999] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_, 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999).
  * Intel [2019] Intel, Math Kernel Library (2019).
  * Wolfram Research, Inc. [2019] Wolfram Research, Inc., Mathematica, Version: 12.0 (2019), Champaign, IL.
  * sci [2019] SciPy, Version: 1.3.0: Open source scientific tools for Python (2019).
  * Hunter [2007] J. D. Hunter, Comput. Sci. Engg. 9, 90 (2007).
  * mat [2019] Matplotlib, Version: 3.1.2 (2019).

## Supplemental Material for
“Cusp singularities in Hessian element distributions of amorphous media”

In this document we provide supplemental figures and derivations related to
the results presented in the main text. In Section 1, we provide detailed
derivations of the partial integrands $\mathcal{P}(H,r)$ that are used in the
computation of the Hessian element distributions $P(H)$ in two and three
dimensions. In Section 2, we display a precise match between the analytic and
numerical distributions obtained using an exponential pair correlation
function $g_{2}(r)$. In Section 3, we provide details of the structure of the
singularities of the partial integrand. In Section 4, we derive exact
asymptotic forms for $P(H)$ in the limit $H\to 0$. Finally in Section 5, we
provide details of the models and software used in our numerical simulations.

###  .1 Hessian Element Distribution

In this Section we derive exact forms for the partial integrands in both two
and three dimensions as announced in Eqs. (6) and (7) in the main text. The
elements of the Hessian matrix may be described in terms of the inter-particle
distance vector ($\mathbf{r}^{ij}$), its magnitude ($r^{ij}$), and the pair-
wise interaction potential ($\psi^{ij}$), which for central potentials is a
function purely of the inter-particle distance $r^{ij}$ (Eq. (3) in main text)

$\mathcal{H}_{\alpha\beta}^{ij}(\mathbf{r}^{ij})=-\left(\frac{\psi_{rr}^{ij}}{{(r^{ij})}^{2}}-\frac{\psi_{r}^{ij}}{{(r^{ij})}^{3}}\right)r^{ij}_{\alpha}r^{ij}_{\beta}-\delta_{\alpha\beta}\frac{\psi_{r}^{ij}}{r^{ij}},$
(13)

where the subscripts of $r$ indicate partial derivatives
$\psi_{r}=\partial_{r^{ij}}\psi$ and
$\psi_{rr}=\partial_{r^{ij}}\partial_{r^{ij}}\psi$. Henceforth, we use $H$ and
$r$ to represent the Hessian element and the inter-particle distance
respectively. The distribution of the elements may then be computed as (Eq.
(4) in main text)

$P\left(H\right)=\int
dr\;d\Omega\;P\left(r,\Omega\right)\delta\left(H-\mathcal{H}_{\alpha\beta}^{ij}(\mathbf{r}^{ij})\right).$
(14)

In addition, assuming a product form for the joint distribution of the inter-
particle distances and angles $P(r,\Omega)=P(r)P(\Omega)$, we have (Eq. (5) in
main text)

$P\left(H\right)=\int_{0}^{r_{c}}dr\;P(r)\;\frac{P\left(\Omega\right)}{\left\lvert\frac{\partial
H}{\partial\Omega}\right\rvert}=\int_{0}^{r_{c}}dr\;P(r)\;\mathcal{P}(H,r).$
(15)

In the sections that follow, we focus our attention on determining the partial
integrand $\mathcal{P}(H,r)$.

#### Diagonal elements

We first consider the diagonal elements $(\alpha=\beta)$ in two dimensions.
Since by symmetry, the $xx$ and $yy$ elements are equivalent, we consider just
the former. The inter-particle distance vector is determined by its magnitude
and angle ($\phi$) with respect to the $x$-axis. Using this we arrive at

$\mathcal{H}_{xx}^{ij}\left(r,\cos\phi\right)=-\left(\psi_{rr}^{ij}-\frac{\psi_{r}^{ij}}{r}\right)\cos^{2}\phi-\frac{\psi_{r}^{ij}}{r}.$
(16)

$-$$\displaystyle\pi$$-$$\displaystyle\pi$/20$\displaystyle\pi$/2$\displaystyle\pi$$\displaystyle\phi$1.01.52.02.53.0$\displaystyle
r$0.00.20.40.60.81.0

Figure 5: Joint distribution of inter-particle distances ($r$) and angles
($\phi$) for energy minimised structures of an R10 system in two dimensions.
This distribution is isotropic in the angular variable $\phi$, justifying a
product form for this joint distribution.

As stated above, we may infer from Fig. 5 that the distributions of $r$ and
$\phi$ are independent, and additionally, that the distribution in $\phi$ is
uniform, with

$\displaystyle P(r,\phi)$ $\displaystyle=P(r)P(\phi),$ $\displaystyle
P_{2\text{D}}(\phi)$ $\displaystyle=\frac{1}{2\pi}.$ (17)

Since the only dependence on $\phi$ is through a cosine, it is convenient to
use the distribution

$P_{2\text{D}}(\cos\phi)=\frac{1}{\pi\sqrt{1-\cos\phi^{2}}}.$ (18)

The partial integrand in this case is

$\mathcal{P}^{2\text{D}}_{\alpha\alpha}={\left[\left\lvert\frac{\partial
H}{\partial\cos\phi}\right\rvert\sqrt{\left\lvert
1-\cos^{2}\phi\right\rvert}\right]}^{-1}.$ (19)

We may use Eq. (16) to simplify this, also noting that

$\cos^{2}\phi=-\frac{H+\frac{\psi_{r}}{r}}{\psi_{rr}-\frac{\psi_{r}}{r}}.$
(20)

Therefore, we have

$\mathcal{P}^{2\text{D}}_{\alpha\alpha}\left(H,r\right)={\left[2\sqrt{\left\lvert\left(H+\frac{\psi_{r}}{r}\right)\left(H+\psi_{rr}\right)\right\rvert}\right]}^{-1},$
(21)

which is the first part of Eq. (6) in the main text.

We next derive the partial integrand in three dimensions. The inter-particle
distance vector is now determined by its magnitude, the polar angle ($\theta$)
subtended on the $z$-axis, and the azimuthal angle ($\phi$) subtended by the
projection of the vector onto the $x-y$ plane, measured relative to the
$x$-axis. By symmetry, the distributions of the $xx$, $yy$ and $zz$ Hessian
elements are the same, and therefore we only need consider one of them

$\mathcal{H}_{zz}^{ij}\left(r^{ij},\theta\right)=-\left(\psi_{rr}^{ij}-\frac{\psi_{r}^{ij}}{r}\right)\cos^{2}\theta-\frac{\psi_{r}^{ij}}{r}.$
(22)

The partial integrand $\mathcal{P}(H,r)$ is then similar to that in two
dimensions (Eq. (19)), since the $zz$ element does not depend on the azimuthal
angle $\phi$. Next, we may infer from Fig. 6 that the distribution of
$\cos\theta$ is uniform

$P_{3\text{D}}(\cos\theta)=\frac{1}{2}.$ (23)

Since the $zz$ element does not involve $\phi$, the integral over $\phi$ in
Eq. (14) evaluates to a constant, and therefore, we arrive at the partial
integrand

$\mathcal{P}^{3\text{D}}_{\alpha\alpha}={\left\lvert\frac{\partial
H}{\partial\cos\theta}\right\rvert}^{-1}.$ (24)

Solving Eq. (22) for $\cos\theta$ as in Eq. (20), and substituting it in Eq.
(24) above, we have

$\mathcal{P}^{3\text{D}}_{\alpha\alpha}\left(H,r\right)={\left[2\sqrt{\left\lvert\left(H+\frac{\psi_{r}}{r}\right)\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\right\rvert}\right]}^{-1},$
(25)

which is the first part of Eq. (7) in the main text.

$-$$\displaystyle\pi$$-$$\displaystyle\pi$/20$\displaystyle\pi$/2$\displaystyle\pi$$\displaystyle\phi$1.21.41.61.82.02.22.4$\displaystyle
r$0.00.20.40.60.8

$-$1.0$-$0.50.00.51.0$\displaystyle\cos(\theta)$1.21.41.61.82.02.22.4$\displaystyle
r$0.00.20.40.60.8

Figure 6: (Left) Joint distribution of inter-particle distances ($r$) and
azimuthal angles ($\phi$) and (Right) joint distribution of inter-particle
distances ($r$) and polar angles ($\cos\theta$) for energy minimised
structures of an R10 system in three dimensions. In both cases the angular
distributions are isotropic.

#### Off-Diagonal elements

The off-diagonal elements ($\alpha\neq\beta$) of the Hessian in two dimensions
may be represented as

$\mathcal{H}_{xy}^{ij}\left(r,\phi\right)=-\left(\psi_{rr}^{ij}-\frac{\psi_{r}^{ij}}{r}\right)\cos\phi\sin\phi.$
(26)

Once again, the partial integrand takes the form

$\mathcal{P}^{2\text{D}}_{\alpha\beta}={\left[\left\lvert\frac{\partial
H}{\partial\cos\phi}\right\rvert\sqrt{\left\lvert
1-\cos^{2}\phi\right\rvert}\right]}^{-1},$ (27)

which can be simplified to yield

$\mathcal{P}^{2\text{D}}_{\alpha\beta}={\left[\left\lvert\frac{\partial
H}{\partial\cos\phi}\sin\phi\right\rvert\right]}^{-1}.$ (28)

Using Eq. (26), the Jacobian factor is given by

$\frac{\partial
H}{\partial\cos\phi}=-\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\sin\phi+\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\frac{\cos^{2}\phi}{\sin\phi}=\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\frac{\cos
2\phi}{\sin\phi}.$ (29)

Substituting this back into Eq. (28), we arrive at

$\mathcal{P}^{2\text{D}}_{\alpha\beta}={\left[\left\lvert\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)\cos
2\phi\right\rvert\right]}^{-1}.$ (30)

We may also rewrite Eq. (26) as

$\sin 2\phi=\frac{-2H}{\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}.$ (31)

Therefore

$\cos
2\phi=\sqrt{1-\frac{4H^{2}}{{\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}^{2}}}.$
(32)

Thus, the partial integrand in terms of $H$ and $r$ simplifies to

$\mathcal{P}^{2\text{D}}_{\alpha\beta}\left(H,r\right)={\left[\sqrt{\left\lvert{\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}^{2}-4H^{2}\right\rvert}\right]}^{-1},$
(33)

which is the second part of Eq. (6) in the main text.

Next, for the calculation for off-diagonal elements in three dimensions, it is
convenient to choose $\mathcal{H}_{xz}$, which now additionally involves the
polar angle $\theta$ as

$\mathcal{H}_{xz}^{ij}\left(r,\theta,\phi\right)=-\left(\psi_{rr}^{ij}-\frac{\psi_{r}^{ij}}{r}\right)\cos\theta\sin\theta\cos\phi.$
(34)

The distribution of $\phi$ in 3D is similar to 2D, as can be seen from Fig. 6,
with $P_{3\text{D}}(\phi)=\frac{1}{2\pi}$. Therefore the partial integrand is
given by

$\mathcal{P}^{3\text{D}}_{\alpha\beta}=\int_{-1}^{1}d(\cos\theta){\left[\left\lvert\frac{\partial
H}{\partial\cos\phi}\right\rvert\sqrt{\left\lvert
1-\cos^{2}\phi\right\rvert}\right]}^{-1}.$ (35)

The Jacobian in the above equation can be evaluated using Eq. (34). Finally,
substituting for $\cos\phi$ using Eq. (34), we have

$\mathcal{P}^{3\text{D}}_{\alpha\beta}=\int_{-1}^{1}d(\cos\theta){\left[\sqrt{\left\lvert{\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}^{2}\sin^{2}\theta\cos^{2}\theta-H^{2}\right\rvert}\right]}^{-1}.$
(36)

We may simplify this integral, and express it as

$I\left(H;\kappa\right)=\frac{\kappa}{H}\int^{1}_{-1}\frac{dx}{\sqrt{x^{2}(1-x^{2})-\kappa^{2}}},$
(37)

which is the second part of Eq. (7) in the main text, and where

$\kappa=\frac{H}{\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}.$ (38)

Since we are interested in the $H\to 0$ limit, we may assume $\kappa<1/2$.
This guarantees that the term within the square root in Eq. (37) remains real.
The limits such that the integral in Eq. (37) is real, are given by

$a=\sqrt{\frac{1-\sqrt{1-4\kappa^{2}}}{2}}\text{ and
}b=\sqrt{\frac{1+\sqrt{1-4\kappa^{2}}}{2}}.$ (39)

We may then write

$I\left(H;\kappa\right)=\frac{\kappa}{H}\int^{b}_{a}\frac{2dx}{\sqrt{x^{2}(1-x^{2})-\kappa^{2}}},$
(40)

Unfortunately, the integral above does not have a closed form solution,
however we can still extract the asymptotic behaviour exactly. In order to
make progress, let us re-write the integral described in Eq. (40), by making
the transformation

$\kappa^{2}=\mu^{2}-\mu^{4},$ (41)

which gives us

$I\left(H;\mu\right)=\frac{\sqrt{\mu^{2}-\mu^{4}}}{H}\int^{\sqrt{1-\mu^{2}}}_{\mu}\frac{2dx}{\sqrt{x^{2}(1-x^{2})-\mu^{2}+\mu^{4}}},$
(42)

where the limits simplify since we may write
$x^{2}\left(1-x^{2}\right)-\mu^{2}+\mu^{4}=\left(x^{2}-\mu^{2}\right)\left(\left(1-x^{2}\right)-\mu^{2}\right)$.

Next, we split the integrand into partial fractions as

$\frac{1}{\sqrt{x^{2}(1-x^{2})-\mu^{2}+\mu^{4}}}=\frac{1}{\sqrt{x^{2}-\mu^{2}}}+\frac{1}{\sqrt{1-x^{2}-\mu^{2}}}+\frac{1-\sqrt{1-x^{2}-\mu^{2}}-\sqrt{x^{2}-\mu^{2}}}{\sqrt{(1-x^{2}-\mu^{2})(x^{2}-\mu^{2})}}.$
(43)

We can now examine the integral over each of these terms separately. For the
first term, we have

$\displaystyle\int^{\sqrt{1-\mu^{2}}}_{\mu}\frac{dx}{\sqrt{x^{2}-\mu^{2}}}$
$\displaystyle\hskip
4.0pt=\left.\tanh^{-1}{\left(\frac{x}{\sqrt{x^{2}-\mu^{2}}}\right)}\right\lvert^{\sqrt{1-\mu^{2}}}_{\mu}$
$\displaystyle\hskip
4.0pt=\tanh^{-1}{\left(\sqrt{\frac{\mu^{2}-1}{2\mu^{2}-1}}\right)}+i\frac{\pi}{2}$
$\displaystyle\overset{\mu\to
0}{\approx}\log{(2)}-\log{(\mu)}+O\left(\mu^{2}\right).$ (44)

Note that this is singular at $\mu=0$. Next, the integral over the second term
yields

$\displaystyle\int^{\sqrt{1-\mu^{2}}}_{\mu}\frac{dx}{\sqrt{1-x^{2}-\mu^{2}}}$
$\displaystyle\hskip
4.0pt=\left.\tan^{-1}{\left(\frac{x}{\sqrt{1-x^{2}-\mu^{2}}}\right)}\right\lvert^{\sqrt{1-\mu^{2}}}_{\mu}$
$\displaystyle\hskip
4.0pt=\frac{\pi}{2}-\tan^{-1}{\left(\frac{\mu}{\sqrt{1-2\mu^{2}}}\right)}$
$\displaystyle\overset{\mu\to
0}{\approx}\frac{\pi}{2}-\mu+O\left(\mu^{3}\right).$ (45)

The integral over the third term cannot be performed exactly, so we first
examine its behaviour as $\mu\to 0$

$\frac{1-\sqrt{1-x^{2}-\mu^{2}}-\sqrt{x^{2}-\mu^{2}}}{\sqrt{(1-x^{2}-\mu^{2})(x^{2}-\mu^{2})}}\overset{\mu\to
0}{\approx}-\frac{1+x-\sqrt{1-x^{2}}}{x(1+x)}+O\left(\mu^{2}\right).$ (46)

We therefore have

$\displaystyle\int^{\sqrt{1-\mu^{2}}}_{\mu}\frac{(-1-x+\sqrt{1-x^{2}})dx}{x(1+x)}$
$\displaystyle\hskip
4.0pt=\sin^{-1}{(\mu)}-\cos^{-1}{(\mu)}-\log\left(\frac{\mu+1}{\sqrt{1-\mu^{2}}+1}\right)$
$\displaystyle\overset{\mu\to 0}{=}\log{(2)}-\frac{\pi}{2}.$ (47)

Putting together the terms from Eqs. (44), (45) and (47), into Eq. (42) and
ignoring terms of $O(\mu)$, we have

$I\left(H;\mu\right)\approx\frac{2\sqrt{\mu^{2}-\mu^{4}}}{H}\left[\log{(4)}-\log{(\mu)}+O\left(\mu\right)\right].$
(48)

Thus the singular behaviour of the full integrand occurs due to the first term
in the partial fractions. Finally, we transform back to $\kappa$ using the
appropriate solution, i.e., the positive real branch in Eq. (41) given by

$\mu=\frac{\sqrt{1-\sqrt{1-4\kappa^{2}}}}{\sqrt{2}}$ (49)

Therefore $I(H;\kappa)$ (Eq. (37)), in the limit of small $\kappa$, simplifies
to

$I\left(H;\kappa\right)\approx\frac{2\kappa}{H}\left[\log{(4)}-\log{(\kappa)}+O\left(\kappa^{2}\right)\right]$
(50)

Thus, the partial integrand for the off-diagonal elements in three dimensions,
when $\kappa\ll 1$ is given by

$\mathcal{P}^{3\text{D}}_{\alpha\beta}\left(H,r\right)=2{\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}^{-1}\log{\left(\frac{4\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}{H}\right)},$
(51)

which is Eq. (8) in the main text.

$\displaystyle{-10^{1}}$$\displaystyle{-10^{-2}}$$\displaystyle{-10^{-5}}$$\displaystyle{-10^{-8}}$$\displaystyle{0}$$\displaystyle{10^{-8}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{1}}$$\displaystyle
H$$\displaystyle{10^{-4}}$$\displaystyle{10^{-3}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{0}}$$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle
P\left(H\right)$Simulation (2D)Analytic
(2D)$\displaystyle\sim{\left|H\right|}^{-1/4}$

$\displaystyle{-10^{1}}$$\displaystyle{-10^{-2}}$$\displaystyle{-10^{-5}}$$\displaystyle{-10^{-8}}$$\displaystyle{0}$$\displaystyle{10^{-8}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{1}}$$\displaystyle
H$$\displaystyle{10^{-4}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle
P\left(H\right)$Simulation (3D)Analytic
(3D)$\displaystyle\sim{\left|H\right|}^{-1/4}$

Figure 7: Distribution of diagonal Hessian elements ($\alpha=\beta$)
corresponding to A-B interactions in the R10 model, using an exact form for
the radial distribution function $g_{2}(r)$ given in Eq. (52) for (Left) two
dimensional systems and (Right) three dimensional systems.

###  .2 $P(H)$ using an exact radial distribution function

The analytical prediction of the Hessian element distribution as shown in Fig.
4 in the main text contains slight discrepancies which occur due to the
interpolated form of the $g_{2}(r)$ used. In this Section we analyse the
distribution of Hessian elements sampled from an artificial, but exact
distribution of pair distances given by

$g_{2}(r)=e^{-(r-1)}\Theta(r-1),$ (52)

where $\Theta(r)$ is the Heaviside theta function. Our results for diagonal
Hessian elements ($\alpha=\beta$) are displayed in Fig. 7 for both two and
three dimensional systems. We find that the analytic procedure yields a
prediction that agrees with the numerical distribution to within numerical
precision.

###  .3 Poles in the complex-$r$ plane

Signs of($\psi_{\delta},H$)$n=3$$n=2$$H^{1/n}$$H^{1/n-1}$(+ ,
+)$Re$$Im$(a)$Re$$Im$(b)(+ , -)$Re$$Im$(c)$Re$$Im$(d) Figure 8: An
illustration of the structure of poles of Eq. (54) in the complex-$(r_{c}-r)$
plane. The figure shows the poles for degrees of smoothness $n=2$ and $n=3$.
The behaviour depends on the relative signs of the potential near the cut-off
$\left(\psi_{\delta}\equiv\psi(r_{c}-\delta)\right)$ and the Hessian element
$H$. The (red) crosses indicate the largest real positive pole that controls
the upper limit of the integral in Eq. (21). The shift $s(H)$ is given by the
distance of this pole to the origin.

In this Section we analyse the variation of the upper limit in $r$ that
appears in the integral expression for $P(H)$ in Eq. (15) (Eq. (5) in the main
text). For example, in the case $H>0$, the contours plotted in the inset of
Fig. 2 in the main text indicate a maximum value of $r$ (i.e., the turn-around
points of the contours). This implies a divergence stemming from the Jacobian
factor in the integral, which can be extracted from the roots of the
denominator of $\mathcal{P}(H,r)$ given in Eqs. (6) and (7). Below, we will
focus on the non-trivial case of diagonal elements in 2D (Eq. (21)). The other
cases are treated in a similar manner, and are significantly simpler.

As stated in Eq. (9) in the main text, for a given degree of smoothness $n$ at
the cut-off $r_{c}$, the interaction potential and its derivatives may be
approximated in terms of the approach to $r_{c}$ as

$\displaystyle\psi\left(r\right)$
$\displaystyle={\left(r_{c}-r\right)}^{n+1}f\left(r\right)$
$\displaystyle\psi_{r}/r$ $\displaystyle\approx
C_{1}{\left(r_{c}-r\right)}^{n}$ $\displaystyle\psi_{rr}$
$\displaystyle\approx C_{2}{\left(r_{c}-r\right)}^{n-1},$ (53)

where $f(r)$ is a regular function, along with $C_{1}=-(n+1)f(r)/r$, and
$C_{2}=n(n+1)f(r)$ which vary significantly slower compared to the power-law
term as $r\to r_{c}$. The singular behaviour in the partial integrand in Eq.
(21) occurs under one of the following conditions:

$\displaystyle{\left(r_{c}-r\right)}^{n}$ $\displaystyle=-\frac{H}{C_{1}}$
$\displaystyle{\left(r_{c}-r\right)}^{n-1}$ $\displaystyle=-\frac{H}{C_{2}}.$
(54)

Using this approximation, for given a value of $H$, the value of $r$ at which
the partial integrand is singular is _shifted_ from the cut-off, which we
quantify as

$r^{*}(H)=r_{c}-s(H).$ (55)

The shifts $s(H)$ depend on the relative signs of the potential near cut-off
$\left(\psi_{\delta}\equiv\psi(r_{c}-\delta)\text{ with }\delta/r_{c}\ll
1\right)$, and the value of the Hessian element $H$. For the sake of
convenience, let us denote the signs of these two quantities as a pair:
$\left(\psi_{\delta},H\right)$.

The value of $r$ at which the partial integrand $\mathcal{P}(H,r)$ becomes
singular is determined by the largest, real, positive pole in Eq. (54). The
shift $s(H)$ can now be identified based on the fact that positive unity
always has a positive root, while negative unity does not. The structure of
the poles of the partial integrand in 2D (for $n=2$ and $3$) are displayed in
Fig. 8, with the positive real roots that contribute to the shift denoted by
(red) crosses. The shifts for these cases have been summarized in Table 2.

Signs of | Shift
---|---
$\left(\psi_{\delta},H\right)$ | $s(H)$
$(+,+)$ | ${(H/C_{1})}^{1/n}$
$(+,-)$ | ${(H/C_{2})}^{1/n-1}$
$(-,+)$ | ${(H/C_{2})}^{1/n-1}$
$(-,-)$ | ${(H/C_{1})}^{1/n}$
Table 2: The various possible shifts $s(H)$ in the upper limit of the integral
in Eq. (15), for the diagonal elements $(\alpha=\beta)$ in two dimensions. The
value of the shift depends on the relative signs of $H\text{ and
}\psi_{\delta}$ as well as the degree of smoothness $n$.

###  .4 Asymptotic Integrals

The behaviour of the partial integrand $\mathcal{P}(H,r)$ as one approaches
the upper limit $r^{*}(H)$, controls the asymptotic form of $P(H)$. It is
therefore convenient to parametrise these expressions in terms of the distance
$\epsilon$ to this singular value $r^{*}(H)$, with

$\epsilon=r^{*}(H)-r=r_{c}-s(H)-r.$ (56)

The partial integrands display different regimes of decay in $\epsilon$ which
are controlled by the distance to the singularities in the complex-$r$ plane
as described in the previous section. For example, in the case of $n=2$, Fig 9
displays the variation of the partial integrand in Eq. (21) at three different
values of $H=10^{-2},10^{-3}$ and $10^{-4}$. In this case, there are two
regimes of decay, with the point of turn-around scaling as $\sqrt{\lvert
H\rvert}$, as may be seen in the scaling collapse displayed in the inset.

$\displaystyle{10^{-6}}$$\displaystyle{10^{-4}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{0}}$$\displaystyle\epsilon$$\displaystyle{10^{-5}}$$\displaystyle{10^{-3}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{1}}$$\displaystyle{10^{3}}$$\displaystyle{10^{5}}$$\displaystyle\mathcal{P}\left(H,\epsilon\right)$$\displaystyle\epsilon^{-1/2}$$\displaystyle\epsilon^{-3/2}$$\displaystyle{10^{-4}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{2}}$$\displaystyle\sqrt{\left|H\right|}\epsilon$$\displaystyle{10^{-8}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{1}}$$\displaystyle{\left|H\right|}^{3/4}\mathcal{P}\left(H,\epsilon\right)$$\displaystyle
H=-10^{-2}$$\displaystyle H=-10^{-4}$$\displaystyle H=-10^{-6}$

Figure 9: Plot of the asymptotic partial integrand $\mathcal{P}(H,\epsilon)$
given in Eq. (59), corresponding to A-B interactions of the R10 model with
degree of smoothness $n=2$. The two regimes of asymptotic behaviour are given
by $[0,\sqrt{H})$ and $[\sqrt{H},\infty)$, with decays $\epsilon^{-1/2}$ and
$\epsilon^{-3/2}$ respectively, as predicted in Eqs. (61) and (60). (Inset)
Scaling collapse of the integrand with the scaled variable $\sqrt{\lvert
H\rvert}\epsilon$.

We now examine each of the partial integrands derived in Sec. .1, in the limit
of the inter-particle distances approaching the singular value $r^{*}(H)$,
i.e. $\epsilon\ll 1$, as described in Sec. .3. We first derive limiting forms
for the partial integrands, identifying the various regimes of decay in
$\epsilon$. Finally we use these asymptotic forms to perform the integration
in Eq. (15), to extract the dominant behaviour of $P(H)$ as $H\to 0$. In terms
of the distance $\epsilon$ and the shift $s\equiv s(H)$, we may rewrite Eq.
(53) as

$\psi_{r}/r\approx C_{1}{\left(\epsilon+s\right)}^{n},\quad\psi_{rr}\approx
C_{2}{\left(\epsilon+s\right)}^{n-1},$ (57)

which is Eq. (11) in the main text. In the rest of this section, we focus on
the case when the potential is positive near cut-off $(\psi_{\delta}>0)$, as
is the case for the R10 model that we study numerically (displayed in Fig. 2
in the main text). The case for $(\psi_{\delta}<0)$ proceeds exactly in the
same manner, with $H\to-H$. Note that the signs of the derivatives of the
potential near the cut-off, alternate with the order of derivative. This may
be seen from Eq. (53), and implies that $\psi_{\delta},C_{1}\text{ and }C_{2}$
have alternately opposing signs. This allows us to determine the shifts from
Eq. (54), purely in terms of the signs of $H$, given the sign of
$\psi_{\delta}$.

#### Diagonal elements

Two Dimensions:  For the case of diagonal elements in two dimensions,
inserting Eq. (57) in Eq. (21), we have

$\mathcal{P}^{2\text{D}}_{\alpha\alpha}\left(H,\epsilon\right)\sim{\left\\{\left[H+C_{1}{(\epsilon+s)}^{n}\right]\left[H+C_{2}{(\epsilon+s)}^{n-1}\right]\right\\}}^{-1/2}.$
(58)

As detailed in Section .3, the value of the shift $s(H)$ depends on the sign
of $H$, i.e., the direction of approach to $0$.

(a) Let us first consider the case of $H<0$. As shown in Fig. 8 (c,d) and
Table 2, for this case, $s\sim{\lvert H\rvert}^{1/(n-1)}$. Eq. (58), ignoring
an overall scale factor, may then be written as

$\mathcal{P}^{2\text{D}}_{\alpha\alpha}\left(H<0,\epsilon\right)\sim{\left\\{\left[-\lvert
H\rvert-{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n}\right]\left[-\lvert
H\rvert+{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n-1}\right]\right\\}}^{-\frac{1}{2}}.$ (59)

At small values of $\epsilon$, this asymptotic form simplifies into power-laws
in $\epsilon$. The number of power-law regimes present are governed by the
distance to the poles in the complex-$\epsilon$ plane. As an illustrative
example, let us focus on the case of $n=2$. The structure of the poles for
this case have been illustrated in Fig. 8 (d). Here, all the poles are at the
same distance scale $\sqrt{\lvert H\rvert}$ away from $\epsilon=0$, implying
two regimes of decay in $\epsilon$, separated by this scale. We illustrate
this in Fig. 9 where we plot the asymptotic partial integrand in Eq. (59). The
inset shows a scaling collapse of the partial integrand for different values
of $H$, displaying the $\sqrt{\lvert H\rvert}$ scale in the turning point.

Finally, substituting Eq. (59) into the integral for $P(H)$ (Eq. (15)), we can
identify the contributions of each of the slopes to the asymptotic behaviour
of $P(H)$ as $H\to 0$. Below, we present these simplifications and
integrations for each regime, with the dominant terms marked with boxes. For
the case of $n=2$, we split the integral in Eq. (59) into two parts $I_{1}(H)$
and $I_{2}(H)$, arising from the regions $\epsilon\in[\sqrt{H},\infty)$ and
$\epsilon\in[0,\sqrt{H})$ respectively. In the first case, the singular
contribution to $P(H)$ is given by

$I_{1}\left(H\right)\sim\int_{\sqrt{H}}^{\infty}{\left[(\epsilon^{2})(\epsilon)\right]}^{-\frac{1}{2}}=\int_{\sqrt{H}}^{\infty}\epsilon^{-\frac{3}{2}}=\left.\epsilon^{-\frac{1}{2}}\right\rvert^{\infty}_{\sqrt{H}}\sim\boxed{H^{-\frac{1}{4}}}.$
(60)

Similarly, when $\epsilon\in[0,\sqrt{H})$, the singular contribution is given
by

$I_{2}\left(H\right)\sim\int_{0}^{\sqrt{H}}{\left[(H)(\epsilon)\right]}^{-\frac{1}{2}}=H^{-\frac{1}{2}}\int_{0}^{\sqrt{H}}\epsilon^{-\frac{1}{2}}=H^{-\frac{1}{2}}\left.\epsilon^{\frac{1}{2}}\right\rvert_{0}^{\sqrt{H}}=H^{-\frac{1}{4}}.$
(61)

Note that the analysis above, predicts power laws of $\epsilon^{-3/2}$ and
$\epsilon^{-1/2}$ of the integrand, as displayed in Fig. 9. Thus, both regimes
contribute to the dominant singular behaviour in $P(H)$ which has the
singularity $H^{-\frac{1}{4}}$, as displayed in Fig. 10. Next, for $n\geq 3$,
an additional (smaller) distance scale also appears, namely $(\sim{\lvert
H\rvert}^{1/n-1})$, as shown in Fig. 8 (c). This implies three regimes of
decay in $\epsilon$, separated by the two scales ${\lvert H\rvert}^{1/n-1}$
and ${\lvert H\rvert}^{1/n}$. Once again we can evaluate the contribution from
each of these regimes to the singular behaviour in $P(H)$ and identify the
dominant singularity.

We treat the case $n=3$ separately, as this presents logarithmic singularities
as we show below. In this case we split the integral in Eq. (59) into three
parts

$I_{1}\left(H\right)\sim\int_{{H}^{\frac{1}{3}}}^{\infty}{\left[(\epsilon^{3})(\epsilon^{2})\right]}^{-\frac{1}{2}}=\int_{{H}^{\frac{1}{3}}}^{\infty}\epsilon^{-\frac{5}{2}}=\left.\epsilon^{-\frac{3}{2}}\right\rvert_{{H}^{\frac{1}{3}}}^{\infty}=H^{-\frac{1}{2}}.$
(62)
$I_{2}\left(H\right)\sim\int_{{H}^{\frac{1}{2}}}^{{H}^{\frac{1}{3}}}{\left[(H)(\epsilon^{2})\right]}^{-\frac{1}{2}}=H^{-\frac{1}{2}}\int_{{H}^{\frac{1}{2}}}^{{H}^{\frac{1}{3}}}\epsilon^{-1}=H^{-\frac{1}{2}}\left.{\log{\epsilon}}\right\rvert_{{H}^{\frac{1}{2}}}^{{H}^{\frac{1}{3}}}\sim\boxed{H^{-\frac{1}{2}}\log{H^{-1}}}.$
(63) $I_{3}\left(H\right)\sim\int_{0}^{\sqrt{H}}{\left[(H)(\epsilon
H^{\frac{1}{2}})\right]}^{-\frac{1}{2}}=H^{-\frac{3}{4}}\int_{0}^{\sqrt{H}}\epsilon^{-\frac{1}{2}}=H^{-\frac{3}{4}}\left.\epsilon^{\frac{1}{2}}\right\rvert_{0}^{\sqrt{H}}=H^{-\frac{1}{2}}.$
(64)

Therefore the dominant singular behaviour in $P(H)$ for $n=3$ is accompanied
by a logarithmic singularity. This behaviour has been displayed using results
from simulations in Fig. 10. Next, we derive results for any general $n>3$.
Once again, splitting the integral in Eq. (59) into three parts, we have

$I_{1}\left(H\right)\sim\int_{{H}^{\frac{1}{n}}}^{\infty}{\left[(\epsilon^{n})(\epsilon^{n-1})\right]}^{-\frac{1}{2}}=\int_{{H}^{\frac{1}{n}}}^{\infty}\epsilon^{-n+\frac{1}{2}}=\left.\epsilon^{-n+\frac{3}{2}}\right\rvert_{{H}^{\frac{1}{n}}}^{\infty}=H^{-1+\frac{3}{2n}}.$
(65)
$I_{2}\left(H\right)\sim\int_{{H}^{\frac{1}{n-1}}}^{{H}^{\frac{1}{n}}}{\left[(H)(\epsilon^{n-1})\right]}^{-\frac{1}{2}}=H^{-\frac{1}{2}}\int_{{H}^{\frac{1}{n-1}}}^{{H}^{\frac{1}{n}}}\epsilon^{-\frac{n}{2}+\frac{1}{2}}=H^{-\frac{1}{2}}\left.{\epsilon^{-\frac{n}{2}+\frac{3}{2}}}\right\rvert_{{H}^{\frac{1}{n-1}}}^{{H}^{\frac{1}{n}}}\sim\boxed{H^{-1+\frac{1}{n-1}}}.$
(66) $I_{3}\left(H\right)\sim\int_{0}^{{H}^{\frac{1}{n-1}}}{\left[(H)(\epsilon
H^{\frac{n-2}{n-1}})\right]}^{-\frac{1}{2}}=H^{-\frac{2n-3}{2n-2}}\int_{0}^{{H}^{\frac{1}{n-1}}}\epsilon^{-1/2}=H^{-\frac{2n-3}{2n-2}}\left.\epsilon^{1/2}\right\rvert_{0}^{{H}^{\frac{1}{n-1}}}\sim
H^{-1+\frac{1}{n-1}}.$ (67)

Above we have simplified the expressions in $I_{3}$ via a binomial expansion
in the second term of Eq. (59). In this case, we find two competing powers
arising from the three regions. The boxed result is dominant for $n>3$,
whereas the contribution of the first region is dominant in the $n\leq 3$
cases as we have detailed above. We verify this non-trivial _switching_
behaviour with increasing degree of smoothness $n$ with results from
simulations that are displayed in Fig. 10.

$\displaystyle{-10^{1}}$$\displaystyle{-10^{-2}}$$\displaystyle{-10^{-5}}$$\displaystyle{-10^{-8}}$$\displaystyle{0}$$\displaystyle{10^{-8}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{1}}$$\displaystyle
H$$\displaystyle{10^{-8}}$$\displaystyle{10^{-6}}$$\displaystyle{10^{-4}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle{10^{4}}$$\displaystyle
P^{2D}_{\alpha\alpha}\left(H\right)$$\displaystyle n=$ 2$\displaystyle n=$
3$\displaystyle n=$
4$\displaystyle\sim{\left|H\right|}^{-1/4}$$\displaystyle\sim{\left|H\right|}^{-1/2}$$\displaystyle\sim{\left|H\right|}^{-1/2}$$\displaystyle\log{\left|H\right|}^{-1}$$\displaystyle\sim{\left|H\right|}^{-5/8}$$\displaystyle\sim{\left|H\right|}^{-2/3}$

Figure 10: Numerically sampled distribution (points) of diagonal Hessian
elements $(\alpha=\beta)$ in energy minimized configurations of the R10 system
in two dimensions, for various degrees of smoothness ($n$) of the potentials
at cut-off. The plots show Hessian elements corresponding only to A-B
interactions. The asymptotic behaviour of the singularities match our
theoretical predictions (lines) given in Eqs. (60), (63) and (66).

(b) Next, we consider the case of $H>0$. In this case, the pole structure
displayed in Fig. 8 (a,b) and detailed in Table 2, yields a shift
$s\sim{\lvert H\rvert}^{1/n}$. Once again, Eq. (58), up to an overall scale
factor, can be simplified to yield

$\mathcal{P}^{2\text{D}}_{\alpha\alpha}\left(H>0,\epsilon\right)\sim{\left\\{\left[\lvert
H\rvert-{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n}}\right)}^{n}\right]\left[\lvert
H\rvert+{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n-1}\right]\right\\}}^{-\frac{1}{2}}.$ (68)

In this case, there is only one distance scale away from the real positive
pole: $\sim{\lvert H\rvert}^{1/n}$, and consequently only two regimes of decay
in $\epsilon$. Following the same procedure as for the case of $H<0$, we split
the integral in Eq. (15) into two parts, yielding the asymptotic behaviours

$I_{1}\left(H\right)\sim\int_{{H}^{\frac{1}{n}}}^{\infty}{\left[(\epsilon^{n})(\epsilon^{n-1})\right]}^{-\frac{1}{2}}=\int_{{H}^{\frac{1}{n}}}^{\infty}\epsilon^{-n+\frac{1}{2}}=\left.\epsilon^{-n+\frac{3}{2}}\right\rvert_{{H}^{\frac{1}{n}}}^{\infty}=\boxed{H^{-1+\frac{3}{2n}}}.$
(69) $I_{2}\left(H\right)\sim\int_{0}^{{H}^{\frac{1}{n}}}{\left[(\epsilon
H^{\frac{n-1}{n}})(H^{\frac{n-1}{n}})\right]}^{-\frac{1}{2}}=H^{-1+\frac{1}{2n}}\int_{0}^{{H}^{\frac{1}{n}}}\epsilon^{-1/2}=H^{-1+\frac{1}{2n}}\left.\epsilon^{1/2}\right\rvert_{0}^{{H}^{\frac{1}{n}}}\sim
H^{-1+\frac{3}{2n}}.$ (70)

Above, the expression in $I_{2}$ has been simplified using a binomial
expansion of the first term in Eq. (68). In this case we find that both
regions contribute the same power $H^{-1+\frac{3}{2n}}$, as corroborated with
simulation results displayed in Fig. 10. A striking aspect of our analysis, is
the prediction and numerical confirmation of the differences in powers of the
singularities, on the negative and positive branches of $P(H)$ for the cases
$n\geq 3$. On the positive side $H>0$ of $P(H)$, there are no changes in the
slope scenarios, with the degree of smoothness $n$.

Three Dimensions:  We next consider the distribution of diagonal Hessian
elements in three dimensions. Inserting Eq. (57) into the expression for the
partial integrand in Eq. (25), we arrive at

$\mathcal{P}^{3\text{D}}_{\alpha\alpha}\left(H,\epsilon\right)\sim{\left\\{\left[H+C_{1}{(\epsilon+s)}^{n}\right]\left[C_{2}{(\epsilon+s)}^{n-1}-C_{1}{(\epsilon+s)}^{n}\right]\right\\}}^{-1/2}.$
(71)

$\displaystyle{-10^{1}}$$\displaystyle{-10^{-2}}$$\displaystyle{-10^{-5}}$$\displaystyle{-10^{-8}}$$\displaystyle{0}$$\displaystyle{10^{-8}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{1}}$$\displaystyle
H$$\displaystyle{10^{-6}}$$\displaystyle{10^{-4}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle{10^{4}}$$\displaystyle
P^{3D}_{\alpha\alpha}\left(H\right)$$\displaystyle n=$ 2$\displaystyle n=$
3$\displaystyle n=$
4$\displaystyle\sim{\left|H\right|}^{-1/4}$$\displaystyle\sim{\left|H\right|}^{-1/2}$$\displaystyle\sim{\left|H\right|}^{-1/2}$$\displaystyle\log{\left|H\right|}^{-1}$$\displaystyle\sim{\left|H\right|}^{-5/8}$$\displaystyle\sim{\left|H\right|}^{-2/3}$

Figure 11: Numerically sampled distribution (points) of diagonal Hessian
elements $(\alpha=\beta)$ in energy minimized configurations of the R10 system
in three dimensions, for various degrees of smoothness ($n$) of the potentials
at cut-off. The plots show Hessian elements corresponding only to A-B
interactions. The asymptotic behaviour of the singularities match our
theoretical predictions (lines) obtained from Eqs. (73), (75) and (76).

(a) As in two-dimensions, let us begin by considering the case of $H<0$.
However, unlike in two dimensions, we now have only one ‘ring’ of poles, since
only the first term in Eq. (71) contains an $H$. In this case, there are no
poles on the positive real axis, implying a shift $s=0$. However, an
‘effective’ shift in the limits of the integral in Eq. (15) still occurs in
order to maintain a real valued integrand. The limits of the integral in $r$
are thus determined by imposing the condition $\cos^{2}\theta\leq 1$ in Eq.
(22). This may also be inferred from the $H<0$ contours in the inset of Fig. 2
in the main text. Using the asymptotic approximations (Eq. (53)) in the afore
mentioned condition yields a shift in the integration limit of
$s(H)\sim{\lvert H\rvert}^{1/(n-1)}$. We therefore arrive at

$\mathcal{P}^{3\text{D}}_{\alpha\alpha}\left(H<0,\epsilon\right)\sim{\left\\{\left[-\lvert
H\rvert-{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n}\right]\left[{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n-1}+{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n}\right]\right\\}}^{-\frac{1}{2}}.$ (72)

Next, let us examine the different regimes of behaviour of the partial
integrand in $\epsilon$. The first region in this case is $\epsilon\in[{\lvert
H\rvert}^{1/n},\infty)$, yielding a contribution

$I_{1}\left(H\right)\sim\int_{{H}^{\frac{1}{n}}}^{\infty}{\left[(\epsilon^{n})(\epsilon^{n-1})\right]}^{-\frac{1}{2}}=\int_{{H}^{\frac{1}{n}}}^{\infty}\epsilon^{-n+\frac{1}{2}}=\left.\epsilon^{-n+\frac{3}{2}}\right\rvert_{{H}^{\frac{1}{n}}}^{\infty}=\boxed{H^{-1+\frac{3}{2n}}}.$
(73)

The second region is $\epsilon\in[{\lvert H\rvert}^{1/(n-1)},{\lvert
H\rvert}^{1/n})$, which yields a contribution (for $n\neq 3$)

$I_{2}\left(H\right)\sim\int_{{H}^{\frac{1}{n-1}}}^{{H}^{\frac{1}{n}}}{\left[(H)(\epsilon^{n-1})\right]}^{-\frac{1}{2}}=H^{-\frac{1}{2}}\int_{{H}^{\frac{1}{n-1}}}^{{H}^{\frac{1}{n}}}\epsilon^{-\frac{n}{2}+\frac{1}{2}}=H^{-\frac{1}{2}}\left.{\epsilon^{-\frac{n}{2}+\frac{3}{2}}}\right\rvert_{{H}^{\frac{1}{n-1}}}^{{H}^{\frac{1}{n}}}=H^{-\frac{1}{2}}\left(H^{-\frac{n-3}{2n}}-H^{-\frac{n-3}{2(n-1)}}\right).$
(74)

Note that for the case $n=3$, the above integral yields a logarithmic
singularity. The dominant term in the integral therefore depends on the value
of the smoothness $n$ as

$\displaystyle I_{2}\left(H\right)\sim\begin{cases}{\left\lvert
H\right\rvert}^{-1+\frac{3}{2n}}&n\in\text{$[2,3)$}\\\ \boxed{{\left\lvert
H\right\rvert}^{\frac{1}{2}}\log{H^{-1}}}&n=3\\\ {\left\lvert
H\right\rvert}^{-1+\frac{1}{n-1}}&n\in(3,\infty).\end{cases}$ (75)

The third region is given by $\epsilon\in[0,{\lvert H\rvert}^{1/(n-1)})$, with
a contribution

$I_{3}\left(H\right)\sim\int_{0}^{{H}^{\frac{1}{n-1}}}{\left[(H)(H)\right]}^{-\frac{1}{2}}=H^{-1}\int_{0}^{{H}^{\frac{1}{n-1}}}d\epsilon=H^{-1}\left.\epsilon\right\rvert_{0}^{{H}^{\frac{1}{n-1}}}\sim\boxed{H^{-1+\frac{1}{n-1}}}.$
(76)

Once again, the behaviour of the dominant slope follows that of $I_{2}$ (Eq.
(75)). This behaviour is illustrated in the results from simulations shown in
Fig. 11.

(b) We next turn to the positive branch of the distribution diagonal elements
in 3D when $H>0$. Similar to the case of two dimensions, the shift is
calculated from the first term in Eq. (71) and is given by $s\sim{\lvert
H\rvert}^{1/n}$, allowing the expression to be simplified to

$\mathcal{P}^{3\text{D}}_{\alpha\alpha}\left(H>0,\epsilon\right)\sim{\left\\{\left[\lvert
H\rvert-{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n}\right]\left[{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n-1}+{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n}\right]\right\\}}^{-\frac{1}{2}}.$ (77)

Similar to the corresponding case in two dimensions, there are no changes in
the slopes with the smoothness. The distance scale is now given by
$\sim{\lvert H\rvert}^{1/n}$. Splitting the integral in Eq. (15) into two
regions, we have the contributions

$I_{1}\left(H\right)\sim\int_{{H}^{\frac{1}{n}}}^{\infty}{\left[(\epsilon^{n})(\epsilon^{n-1})\right]}^{-\frac{1}{2}}=\int_{{H}^{\frac{1}{n}}}^{\infty}\epsilon^{-n+\frac{1}{2}}=\left.\epsilon^{-n+\frac{3}{2}}\right\rvert_{{H}^{\frac{1}{n}}}^{\infty}=\boxed{H^{-1+\frac{3}{2n}}}.$
(78) $I_{2}\left(H\right)\sim\int_{0}^{{H}^{\frac{1}{n}}}{\left[(\epsilon
H^{\frac{n-1}{n}})(H^{\frac{n-1}{n}})\right]}^{-\frac{1}{2}}=H^{-1+\frac{1}{2n}}\int_{0}^{{H}^{\frac{1}{n}}}\epsilon^{-1/2}=H^{-1+\frac{1}{2n}}\left.\epsilon^{1/2}\right\rvert_{0}^{{H}^{\frac{1}{n}}}\sim
H^{-1+\frac{3}{2n}}.$ (79)

The first term in $I_{2}$, has been simplified via a binomial expansion of the
second term in Eq. (77). In this case, both slopes contribute the same power
$H^{-1+\frac{3}{2n}}$,as shown in results from simulations in Fig. 11. As in
two dimensions, there is a difference in slopes on the negative and positive
branches of $P(H)$ in three dimensions.

$\displaystyle{-10^{1}}$$\displaystyle{-10^{-2}}$$\displaystyle{-10^{-5}}$$\displaystyle{-10^{-8}}$$\displaystyle{0}$$\displaystyle{10^{-8}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{1}}$$\displaystyle
H$$\displaystyle{10^{-9}}$$\displaystyle{10^{-7}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-3}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{1}}$$\displaystyle{10^{3}}$$\displaystyle
P^{2D}_{\alpha\beta}\left(H\right)$$\displaystyle n=$ 2$\displaystyle n=$
3$\displaystyle n=$
4$\displaystyle\sim{\left|H\right|}^{-1/2}$$\displaystyle\sim{\left|H\right|}^{-2/3}$$\displaystyle\sim\log{\left|H\right|}^{-1}$$\displaystyle{-10^{1}}$$\displaystyle{-10^{-3}}$$\displaystyle{0}$$\displaystyle{10^{-3}}$$\displaystyle{10^{1}}$0.0000.0250.0500.0750.100

Figure 12: Numerically sampled distribution (points) of off-diagonal Hessian
elements $(\alpha\neq\beta)$ in energy minimized configurations of the R10
system in two dimensions, for various degrees of smoothness ($n$) of the
potentials at cut-off. The plots show Hessian elements corresponding only to
A-B interactions. The asymptotic behaviour of the singularities match our
theoretical predictions (lines) obtained from Eq. (83). (Inset) Semi-log plot
of the distribution for $n=2$, highlighting the logarithmic singularity.

#### Off-Diagonal elements

Two Dimensions:  Next, we turn to the distribution of off-diagonal Hessian
elements in two dimensions, with the partial integrand given in Eq. (33). It
is clear from this expression that $H<0$ and $H>0$, have exactly the same
behaviour as the element $H$ only appears as a quadratic form. Using Eq. (57)
in the partial integrand for off-diagonal elements (Eq. (33)), we have

$\mathcal{P}^{2\text{D}}_{\alpha\beta}\left(H,\epsilon\right)\sim{\left\\{{\left[C_{2}{(\epsilon+s)}^{n-1}-C_{1}{(\epsilon+s)}^{n}\right]}^{2}-4H^{2}\right\\}}^{-1/2}.$
(80)

Computing the relevant pole in the above expression, we find that the shift
$s\sim{\lvert H\rvert}^{1/(n-1)}$. Therefore, the partial integrand, up to an
overall scale factor is given by

$\mathcal{P}^{2\text{D}}_{\alpha\beta}\left(H,\epsilon\right)\sim{\left\\{{\left[{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n-1}+{\left(\epsilon+{\lvert
H\rvert}^{\frac{1}{n-1}}\right)}^{n}\right]}^{2}-H^{2}\right\\}}^{-\frac{1}{2}}.$
(81)

Once again we compute the powers of the singularity in $P(H)$. In this case
the partial integrand has three regimes of decay, that can be extracted by
appropriate binomial expansions of Eq. (81), yielding the regions (i)
$[{\lvert H\rvert}^{1/(n-1)},\infty)$, (ii) $[{\lvert
H\rvert}^{2/(n-1)},{\lvert H\rvert}^{1/(n-1)})$ and (iii) $[0,{\lvert
H\rvert}^{2/(n-1)})$ In region (i), we have

$I_{1}\left(H\right)\sim\int_{{H}^{\frac{1}{n-1}}}^{\infty}{\left[{(\epsilon^{n-1})}^{2}\right]}^{-\frac{1}{2}}=\int_{{H}^{\frac{1}{n-1}}}^{\infty}\epsilon^{-n+1}=\left.\epsilon^{-n+2}\right\rvert_{{H}^{\frac{1}{n-1}}}^{\infty},$
(82)

yielding the following singular behaviours

$\displaystyle I_{1}\left(H\right)\sim\begin{cases}\boxed{\log{{\left\lvert
H\right\rvert}^{-1}}}&n=2\\\ \boxed{{\left\lvert
H\right\rvert}^{-1+\frac{1}{n-1}}}&n\in(2,\infty).\end{cases}$ (83)

The integral in the second region $\epsilon\in[{\lvert
H\rvert}^{2/(n-1)},{\lvert H\rvert}^{1/(n-1)})$ can be performed by binomially
expanding the terms of $(\epsilon+s)$. Therefore the contribution from this
region is

$I_{2}\left(H\right)\sim\int_{{H}^{\frac{2}{n-1}}}^{{H}^{\frac{1}{n-1}}}{\left[{\left(H+\epsilon
H^{\frac{n-2}{n-1}}\right)+\left(H^{\frac{n}{n-1}}\right)}^{2}-H^{2}\right]}^{-\frac{1}{2}}\sim
H^{\frac{-(2n-3)}{2(n-1)}}\int_{{H}^{\frac{2}{n-1}}}^{{H}^{\frac{1}{n-1}}}\epsilon^{-\frac{1}{2}}\sim
H^{-1+\frac{1}{n-1}}$ (84)

The contribution from the third region $\epsilon\in[0,{\lvert
H\rvert}^{2/(n-1)})$ can be evaluated as

$I_{3}\left(H\right)\sim\int_{0}^{{H}^{\frac{2}{n-1}}}{\left[{\left(H+H^{\frac{n}{n-1}}\right)}^{2}-H^{2}\right]}^{-\frac{1}{2}}\sim
H^{\frac{-(2n-1)}{2(n-1)}}\int_{0}^{{H}^{\frac{2}{n-1}}}d\epsilon\sim
H^{-1+\frac{3}{2(n-1)}}$ (85)

Therefore the dominant contribution arises from the first two regions
$(I_{1},I_{2})$. Additionally our analysis predicts, a logarithmic singularity
for the case $n=2$. These behaviours are illustrated in our numerically
sampled distributions of off-diagonal elements in two dimensions in Fig. 12.

$\displaystyle{-10^{1}}$$\displaystyle{-10^{-2}}$$\displaystyle{-10^{-5}}$$\displaystyle{-10^{-8}}$$\displaystyle{0}$$\displaystyle{10^{-8}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{1}}$$\displaystyle
H$$\displaystyle{10^{-6}}$$\displaystyle{10^{-4}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle{10^{4}}$$\displaystyle
P^{3D}_{\alpha\beta}\left(H\right)$$\displaystyle n=$ 2$\displaystyle n=$
3$\displaystyle n=$
4$\displaystyle\sim{\left|H\right|}^{-1/2}$$\displaystyle\sim{\left|H\right|}^{-2/3}$$\displaystyle\sim\log{\left|H\right|}^{-1}$$\displaystyle{-10^{0}}$$\displaystyle{-10^{-4}}$$\displaystyle{0}$$\displaystyle{10^{-4}}$$\displaystyle{10^{0}}$0.00.51.01.5

Figure 13: Numerically sampled distribution (points) of off-diagonal Hessian
elements $(\alpha\neq\beta)$ in energy minimized configurations of the R10
system in three dimensions, for various degrees of smoothness of the
potentials at cut-off. The plots show Hessian elements corresponding only to
A-B interactions. The asymptotic behaviour of the singularities match our
theoretical predictions (lines) obtained from Eq. (88). Semi-log plot of the
distribution for $n=2$, highlighting the logarithmic singularity.

Three Dimensions:  Finally, to find the asymptotic behaviour of off-diagonal
elements in three dimensions, we insert Eq. (57) into the expression for the
partial integrand in the limit
$\kappa=H{\left(\psi_{rr}-\frac{\psi_{r}}{r}\right)}^{-1}\to 0$ in Eq. (51).
We therefore arrive at (ignoring an overall multiplicative constant)

$\mathcal{P}^{3\text{D}}_{\alpha\beta}\left(H,\epsilon\right)\sim{\left(\epsilon^{n-1}+\epsilon^{n}\right)}^{-1}\log{\left\lvert\frac{\epsilon^{n-1}-\epsilon^{n}}{H}\right\rvert}.$
(86)

In this case, the limits of the integral in Eq. (15) are determined by the
region of validity $(\kappa<1)$ of the partial integrand Eq. (37). Since we
are working with the asymptotic expression Eq. (51), there is only one region
to consider, and we have the following contribution to the Hessian element
distribution

$\displaystyle I\left(H\right)$
$\displaystyle\sim\int_{{H}^{\frac{1}{n-1}}}^{\infty}\epsilon^{-n+1}\log{\left(H^{-1}\epsilon^{n-1}\right)}$
$\displaystyle\sim\left.{\frac{\epsilon^{-n+2}}{{(n-2)}^{2}}\left[(n-1)-(n-2)\log{\left(H^{-1}\epsilon^{n-1}\right)}\right]}\right|^{\infty}_{H^{\frac{1}{n-1}}}$
(87)

We therefore obtain the following asymptotic behaviour

$\displaystyle I\left(H\right)\sim\begin{cases}\boxed{\log{{\left\lvert
H\right\rvert}^{-1}}}&n=2\\\ \boxed{{\left\lvert
H\right\rvert}^{-1+\frac{1}{n-1}}}&n\in(2,\infty).\end{cases}$ (88)

These behaviours have been illustrated in the results from our numerical
simulations in Fig. 13. The case of $n=2$ yields a logarithmic singularity, as
may be seen from the inset of the figure.

In summary, we have analysed the singularities associated with the
distributions of Hessian elements for all cases in two and three dimensions,
for any degree of smoothness $n$ of the interaction potential at the cut-off
distance. Remarkably we find the same asymptotic behaviours in the
distribution of Hessian elements for _all_ cases in _both_ two and three
dimensions. This is surprising since the mechanisms leading to these
singularities detailed above are very different. We note that although the
dominant singularity exhibits this universal feature, the sub-dominant
corrections to these behaviours may be different in two and three dimensions.
The power of the singularities of the distributions of diagonal elements
display a change in behaviour at $n=3$, and the off-diagonal elements, at
$n=2$. In both cases, systems at the cross-over value of smoothness, display a
logarithmic singularity in the distribution of Hessian elements.

###  .5 Additional Simulation Details

#### Simulation Potentials

We simulate a 50:50 mixture of two particle types A and B. The interaction
potentials are cut-off at a distance

$r_{c}=1.385418025\;\sigma.$ (89)

with the three interaction diameters given by

$\displaystyle\sigma_{AA}$ $\displaystyle=1.0,$ $\displaystyle\sigma_{BB}$
$\displaystyle=1.4,$ $\displaystyle\sigma_{AB}$
$\displaystyle=\sqrt{\sigma_{AA}\sigma_{BB}}.$ (90)

The only difference between the parameters in two and three dimensions of this
model are the reduced densities given by

$\displaystyle\rho_{2\text{D}}$ $\displaystyle=0.85,$
$\displaystyle\rho_{3\text{D}}$ $\displaystyle=0.81.$ (91)

In our simulations we focus on the purely repulsive pairwise potential, given
by a tenth order polynomial, termed ‘R10’. For the case of degree of
smoothness $n=2$, the potential is given by

$\displaystyle\psi$
$\displaystyle={\left(\frac{\sigma}{r}\right)}^{10}+c_{0}+c_{2}{\left(\frac{r}{\sigma}\right)}^{2}+c_{4}{\left(\frac{r}{\sigma}\right)}^{4},$
$\displaystyle c_{0}=-0.8061409035399235,$ $\displaystyle c_{2}=+0.7,$
$\displaystyle c_{4}=-0.15630021928760743.$ (92)

In the case of the potential smoothed to three derivatives $(n=3)$, we have

$\displaystyle\psi={\left(\frac{\sigma}{r}\right)}^{10}$
$\displaystyle+c_{0}+c_{2}{\left(\frac{r}{\sigma}\right)}^{2}+c_{4}{\left(\frac{r}{\sigma}\right)}^{4}+c_{6}{\left(\frac{r}{\sigma}\right)}^{6},$
$\displaystyle c_{0}$ $\displaystyle=-2.149709076106463,$ $\displaystyle
c_{2}$ $\displaystyle=+2.8000000000000007,$ $\displaystyle c_{4}$
$\displaystyle=-1.2504017543008594,$ $\displaystyle c_{6}$
$\displaystyle=+0.1900092331616435.$ (93)

Finally, in the case of the potential smoothed to four derivatives $(n=4)$, we
have

$\displaystyle\psi={\left(\frac{\sigma}{r}\right)}^{10}+c_{0}$
$\displaystyle+c_{2}{\left(\frac{r}{\sigma}\right)}^{2}+c_{4}{\left(\frac{r}{\sigma}\right)}^{4}+c_{6}{\left(\frac{r}{\sigma}\right)}^{6}+c_{8}{\left(\frac{r}{\sigma}\right)}^{8},$
$\displaystyle c_{0}$ $\displaystyle=-4.836845421239541,$ $\displaystyle
c_{2}$ $\displaystyle=+8.400000000000002,$ $\displaystyle c_{4}$
$\displaystyle=-5.626807894353868,$ $\displaystyle c_{6}$
$\displaystyle=+1.7100830984547917,$ $\displaystyle c_{8}$
$\displaystyle=-0.19798989873223335.$ (94)

#### Software

We use LAMMPS [64, 65] to perform simulations of glass forming liquids and
utilise the built-in conjugate-gradient minimiser with a force tolerance
$\sqrt{\sum_{i=1}^{N}{\lvert F_{i}\rvert}^{2}}\text{ of }10^{-11}$ to reach
inherent structures from an equilibrium sample at a temperature $T=1.0$ (lj
units). Eigenvalue calculations are performed using the dsyevr LAPACK [66]
routine in Intel MKL [67]. Analyses are performed using Mathematica [68] and
SciPy [69]. Plotting is performed using Matplotlib [70, 71].