Cusp singularities in Hessian element distributions of amorphous media

We show that the distribution of elements $H$ in the Hessian matrices associated with amorphous materials exhibit cusp 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 cusp singularities affect the stability of amorphous solids, through the distributions of the minimum eigenvalue of the Hessian matrix.

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][3][4][5][6][7][8][9][10]. An important aspect in the study of such systems is their vibrational properties, typically probed through the eigenvalues of the Hessian matrix [11][12][13], with low-frequency modes being particularly relevant at the low-temperatures where glass physics dominates [13][14][15][16][17][18][19][20][21][22][23][24][25]. Since glass forming systems settle into disordered configurations, their Hessian matrices are naturally characterised by a distribution of elements, which so far has received very little attention [26,27]. Random matrix treatments therefore provide a natural framework with which to model amorphous systems [28][29][30][31][32][33], with many studies having focused on lattices with random interactions [34][35][36]. Even though there are indications of an underlying "amorphous order" in such systems [9,10], their behaviour differs fundamentally from that of crystals, with the low lying eigenvalues of their Hessian matrices displaying marked non-Debye behaviour [21][22][23][24][25], which is related to the growth of a structural length scale [37] as well as the growth of static correlations [38]. Although the Hessian matrix is relatively simple to compute in simulations [39,40], an experimental determination remains difficult, accessible for example, only through displacement correlations in colloidal glasses [41]. In this context, it is important to study the vibrational properties of model systems via simulations.
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 [42]. However, the effect of this smoothness on vibrational properties of such amorphous systems has never been investigated. In crystals, the ordered * vishnuvk@tifrh.res.in † smarajit@tifrh.res.in ‡ kramola@tifrh.res.in 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 cusp 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 smoothnessdependent cusp singularities with a power γ < 0. We derive exact results for γ, for any degree of smoothness, in both two dimensions (2D) as well as three dimensions (3D). Our results are summarized in Table I, 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 {r i }, where the particle index i ∈ {1, . . . , N }. The Hessian is then a dN × dN matrix where d represents the dimension of the system, with elements H ij αβ that describe the stiffness between particles i and j, along the coordinates α, β ∈ {x, y, z}. In terms of the inter-particle distance vector r ij = r i − r j , the Hessian elements may be expressed as Here, U {r i } is the total potential energy of the system, which is a function of the positions of all particles. For pairwise additive interactions, where ψ ij is the interaction potential between particles i and j. We consider central potentials ψ(r), where r ≡ r ij = |r ij | is the distance between particles i and j.
In addition, these potentials are smoothed to n derivatives at a cut-off distance r c , i.e., d m ψ dr m rc = 0 for all 0 ≤ m ≤ 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 [43] H ij where the subscripts of r indicate partial derivatives with respect to the inter-particle distances: Since the potential and its derivatives vanish at the cut-off distance, the small Hessian elements arise primarily due to pairdistances 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 Ω with respect to these fixed axes, which is a function of the angle φ in 2D, and both the polar and azimuthal angles θ, φ in 3D. The distribution of Hessian elements P (H) is then given by where P (r, Ω) 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 The R10 potential ψ(r) ∼ r −10 + c + br 2 + ar 4 for A-B interactions, smoothed to two derivatives (n = 2) at the cut-off (ψ(rc) = ψ (rc) = ψ (rc) = 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 θ), displaying contours which contribute to the distribution of diagonal Hessian elements (α = β) at a fixed H.
function constraint in Eq. (4) selects the corresponding contour in (r, Ω) through Eq. (3), as shown in the inset of Fig. 2. 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, Ω) = P (r)P (Ω). Numerically sampling this joint distribution P (r, Ω) 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 2D (φ) = P 3D (φ) = 1 2π , and P 3D (cos θ) = 1 2 . The distribution of r is the radial distribution function g 2 (r), normalised over the range of interaction r ∈ [0, r c ], i.e., P (r) = g 2 (r)Θ(r c − r)N , where Θ is the Heaviside function, and N −1 = For 3D systems, we choose H zz and H xz to represent the diagonal and off-diagonal elements, respectively. From In the case of diagonal elements, it is convenient to choose Ω ≡ cos θ, with the distribution P 3D (cos θ) = 1 2 (see inset of Fig. 3). Under these conditions, the partial integrand defined in In the case of off-diagonal elements, it is convenient to choose Ω ≡ cos φ, and given that P (φ) = 1 2 , we arrive at the exact integral form Finally, inverting the above expressions for H 3D zz and H 3D xz to express cos θ and cos φ in terms of H and r, we arrive at the simplified expressions with κ = H ψ rr − ψr r −1 . The above integral form for the off-diagonal elements has the asymptotic behaviour (refer to Supplemental Material [44]/ 1) 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 [44]/ 2, we additionally show a precise match between P (H) obtained analytically and numerically, for an exponentially decaying g 2 (r).

Asymptotic Forms and Cusp Singularities:
We focus next on the behaviour of P (H) in the limit H → 0. It is clear from the forms of the partial integrands in Eqs. (6) and (7), that they diverge as r → r c and H → 0, since ψ r → 0 and ψ rr → 0. Consequently, The corresponding integrals determining P (H) in Eq. (5) exhibit singularities P (H) ∼ H γ as H → 0, as displayed in Figs. 1 (b) and 4. We show below that the strength γ 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 → 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 → r c , may therefore be approximated as 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 → r c .
. Under this approximation, the singular points of P(H, r) in the complex-r plane are determined by the two expressions The full structure of the poles in Eq. (10) is detailed in the Supplemental Material [44]/ 3. However, as the  In order to analyse the asymptotic forms of P (H) in Eq. (5) as H → 0, we define a small variable as a distance to r * at a given value of H, = (r c − s) − r. The derivatives of the potential described in Eq. (9) can then be written as 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 ψ δ ≡ ψ(r c − δ) with δ/r c 1. Below, we present the analysis for the diagonal elements in 2D, with the details of all cases presented in the Supplemental Material [44]/ 4. Using Eq. (11) in Eq. (6), for 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 in terms of the behaviour of the partial integrand in Eq.  Table I. Remarkably, although the expressions in Eqs. (6) and (7) have very different forms, they yield exactly the same results for the cusp singularities in both 2D and 3D. 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 [44]/ 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 → 0 + with the exponents − 1 4 for n = 2 and − 1 2 for n = 3, as predicted in Table I. 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 cusp singularity in P (H), by analysing the distribution of the minimum eigenvalue, λ min . Our results 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 cusp singularity that depends on the smoothness of the interaction potential at the cut-off distance. Remarkably, the results for the cusp singularities are exactly the same in both 2D and 3D, a fact that warrants deeper investigation. Our results are particularly relevant for numerical studies of glasses, where the effect of smoothness in interaction potentials has not yet been investigated. 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. Such singularities are also an important ingredient to be taken into account in Random Matrix treatments of Hessian matrices of amorphous materials. It would be interesting to extend our analytic results to construct bounds on the vibrational density of states of amorphous systems [45]. Finally, it would also be interesting to extend our results to study jamming transitions, that are characterized by diverging pair correlation functions at the cut-off distance [46][47][48][49].

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 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 → 0. Finally in Section 5, we provide details of the models and software used in our numerical simulations.

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 (r ij ), its magnitude (r ij ), and the pair-wise interaction potential (ψ ij ), which for central potentials is a function purely of the inter-particle distance r ij (Eq. (3) in main text) where the subscripts of r indicate partial derivatives ψ r = ∂ r ij ψ and ψ rr = ∂ r ij ∂ r ij ψ. 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) In addition, assuming a product form for the joint distribution of the inter-particle distances and angles P (r, Ω) = P (r)P (Ω), we have (Eq. In the sections that follow, we focus our attention on determining the partial integrand P(H, r).

Diagonal elements
We first consider the diagonal elements (α = β) 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 (φ) with respect to the x-axis. Using this we arrive at As stated above, we may infer from Fig. 5 that the distributions of r and φ are independent, and additionally, that the distribution in φ is uniform, with P (r, φ) = P (r)P (φ), Since the only dependence on φ is through a cosine, it is convenient to use the distribution The partial integrand in this case is We may use Eq. (16) to simplify this, also noting that Therefore, we have 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 (θ) subtended on the z-axis, and the azimuthal angle (φ) 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 The partial integrand P(H, r) is then similar to that in two dimensions (Eq. (19)), since the zz element does not depend on the azimuthal angle φ. Next, we may infer from Fig. 6 that the distribution of cos θ is uniform Since the zz element does not involve φ, the integral over φ in Eq. (14) evaluates to a constant, and therefore, we arrive at the partial integrand which is the first part of Eq. (7) in the main text.  6. (Left) Joint distribution of inter-particle distances (r) and azimuthal angles (φ) and (Right) joint distribution of inter-particle distances (r) and polar angles (cos θ) 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 (α = β) of the Hessian in two dimensions may be represented as Once again, the partial integrand takes the form which can be simplified to yield Using Eq. (26), the Jacobian factor is given by Substituting this back into Eq. (28), we arrive at We may also rewrite Eq. (26) as Therefore Thus, the partial integrand in terms of H and r simplifies to 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 H xz , which now additionally involves the polar angle θ as The distribution of φ in 3D is similar to 2D, as can be seen from Fig. 6, with P 3D (φ) = 1 2π . Therefore the partial integrand is given by The Jacobian in the above equation can be evaluated using Eq. (34). Finally, substituting for cos φ using Eq. (34), we have We may simplify this integral, and express it as which is the second part of Eq.
Since we are interested in the H → 0 limit, we may assume κ < 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 We may then write 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 which gives us where the limits simplify since we may write x 2 1 − x 2 − µ 2 + µ 4 = x 2 − µ 2 1 − x 2 − µ 2 . Next, we split the integrand into partial fractions as 1 We can now examine the integral over each of these terms separately. For the first term, we have Note that this is singular at µ = 0. Next, the integral over the second term yields The integral over the third term cannot be performed exactly, so we first examine its behaviour as µ → 0 We therefore have √ Putting together the terms from Eqs. (44), (45) and (47), into Eq. (42) and ignoring terms of O(µ), we have Thus the singular behaviour of the full integrand occurs due to the first term in the partial fractions. Finally, we transform back to κ using the appropriate solution, i.e., the positive real branch in Eq. (41) given by Therefore I(H; κ) (Eq. (37)), in the limit of small κ, simplifies to Thus, the partial integrand for the off-diagonal elements in three dimensions, when κ 1 is given by which is Eq. (8) in the main text.

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 where Θ(r) is the Heaviside theta function. Our results for diagonal Hessian elements (α = β) 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.

Poles in the complex-r plane
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 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 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 → r c . The singular behaviour in the partial integrand in Eq. (21) occurs under one of the following conditions: Signs of (ψ δ , H) n = 3 n = 2 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 The shifts s(H) depend on the relative signs of the potential near cut-off (ψ δ ≡ ψ(r c − δ) with δ/r c 1), 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: (ψ δ , H).
The value of r at which the partial integrand 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 II.

Asymptotic Integrals
The behaviour of the partial integrand 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 to this singular value r * (H), with The partial integrands display different regimes of decay in 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  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.
1, as described in Sec. 3. We first derive limiting forms for the partial integrands, identifying the various regimes of decay in . Finally we use these asymptotic forms to perform the integration in Eq. (15), to extract the dominant behaviour of P (H) as H → 0. In terms of the distance and the shift s ≡ s(H), we may rewrite Eq. (53) as 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 (ψ δ > 0), as is the case for the R10 model that we study numerically (displayed in Fig. 2 in the main text). The case for (ψ δ < 0) proceeds exactly in the same manner, with H → −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 ψ δ , C 1 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 ψ δ .
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 ≤ 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. (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 II, In this case, there is only one distance scale away from the real positive pole: ∼ |H| 1/n , and consequently only two regimes of decay in . 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 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+ 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 cusp singularities, on the negative and positive branches of P (H) for the cases n (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 θ ≤ 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) ∼ |H| 1/(n−1) . We therefore arrive at  Computing the relevant pole in the above expression, we find that the shift s ∼ |H| 1/(n−1) . Therefore, the partial integrand, up to an overall scale factor is given by