Critical energy landscape of linear soft spheres

Disordered systems are often only marginally stable, with an abundance of soft excitations that can push typical configurations on the verge of an instability point. A prominent example of this behavior is provided by random packings of spheres at a special point: jamming. In this work, we show that energy minima of dense soft spheres interacting through a linear ramp potential are critical in an entire phase and display non-linear marginal stability with important features in common with jamming points. Excitations around local minima are non-linear, may be system spanning and are characterized by a set of non-trivial critical exponents. The linear potential makes such configurations isostatic, meaning that they are characterized by an isostatic number of perfectly kissing spheres. As a result, we show that, surprisingly, the critical exponents describing non-linear excitations in the finite energy jammed phase appear to be the same -- within numerical precision -- as the ones appearing at the zero-energy jamming point. The emergence of isostaticity is then accompanied by infinite isostatic lengthscales that characterize the potential energy landscape and make such energy minima critical.

Disordered systems are often only marginally stable, with an abundance of soft excitations that can push typical configurations on the verge of an instability point. A prominent example of this behavior is provided by random packings of spheres at a special point: jamming. In this work, we show that energy minima of dense soft spheres interacting through a linear ramp potential are critical in an entire phase and display non-linear marginal stability with important features in common with jamming points. Excitations around local minima are non-linear, may be system spanning and are characterized by a set of non-trivial critical exponents. The linear potential makes such configurations isostatic, meaning that they are characterized by an isostatic number of perfectly kissing spheres. As a result, we show that, surprisingly, the critical exponents describing non-linear excitations in the finite energy jammed phase appear to be the same -within numerical precisionas the ones appearing at the zero-energy jamming point. The emergence of isostaticity is then accompanied by infinite isostatic lengthscales that characterize the potential energy landscape and make such energy minima critical.
The jamming transition is met when a large collection of hard objects is rapidly compressed until the occupied volume cannot be further decreased. Since more than twenty years, the statistical properties of jamming points have shaped our thinking of low temperature glasses, suggested principles underlying amorphous rigidity, and provided mechanisms to rationalize low energy excitations in glasses [1,2]. As a consequence of the rapid compression, random packings of spheres are just on the verge of mechanical instability: the number of contacts between the spheres is isostatic, meaning that in d dimensions each sphere has on average 2d contacts, that is the least one for which the system can sustain pressure [3][4][5]. While the location of the transition depends on the compression protocol [6], isostaticity is generic. Particles jam when -exactly [7] -one contact per degree of freedom is established [8]. As such, mechanical marginal stability brings about criticality and diverging lengthscales [9]. The jamming point is a critical point characterized by a set of critical exponents describing both the behavior of bulk physical quantities, such as pressure, energy and contacts [4,10] as well as the microstructure of amorphous packings [11][12][13]. In particular, a common characterization of heterogeneity is provided by local statistics of contact forces and interparticle distances. Marginal stability implies power law behavior of the distribution of these quantities at small argument [11,12] and predicts a non trivial relation between the corresponding exponents [14]. These exponents have actually been computed within the exact solution of the glassy phases of spherical particles in the limit of infinite dimensionality [15,16] and have been shown to agree -within numerical precision-with numerical simulations of hard and soft spheres in various physical dimensions [17].
While the jamming point is obviously independent of the interaction potential between spheres, many features of the overcompressed, mechanically stable phase, depend on the nature of the soft potential that is used. If we denote by {x i } i=1,...,N the positions of the centers of N spheres, we can define a gap between two spheres, say i and j, as h ij = r ij − σ ij , where we have denoted by σ ij the sum of the radii of the corresponding spheres and by r ij = |x i − x j | the distance between their centers. In the overcompressed phase, above jamming, one typically defines a pure power . f c is a constant that essentially sets the unit of forces. Common choices for the penalty exponent α are α = 2 or α = 5/2, corresponding respectively to Harmonic and Hertzian spheres. If α > 1 the potential is convex and differentiable at r ij = σ ij , i.e. when spheres just touch. As a consequence, given a contact at jamming, an infinitesimal normal force is enough to destabilize it and can cause an overlap between the corresponding particles. Therefore, for α > 1 jamming is a singular point in the phase diagram: beyond φ J , as soon as the particles overlap, the system stabilizes. It means that jamming criticality is washed out when we enter in the jammed phase, where the intensive energy of the system is finite and the system sits in a local minimum at zero temperature.
In this paper, we investigate the potential energy landscape (PEL) of soft spheres interacting through a linear ramp potential, obtained by setting α = 1, above the jamming transition point. This potential, which is at the boundary between convex and non-convex interparticle potentials, presents important qualitative differences from the case α > 1. First of all, it is non-differentiable: small forces applied to contacts do not necessarily destabilize them. To induce an overlap, a total force greater than f c is necessary. In addition, the modulus of the force generated by an overlap does not depend on the extent of the overlap itself.
We focus on systems of two and three dimensional polydisperse spheres and produce local minima by gradient descent minimization. Our main result is that jamming criticality survives even when we enter in the overcompressed phase where the system sits in a finite energy local minimum of the PEL.
Our main findings are that • accessible local minima of the PEL are isostatic. Even if there is a finite fraction of overlapping spheres making the total energy finite, there is also an isostatic number of pairs of spheres that just touch. We call interacting spheres couples of spheres that either are in perfect contact (contacts) h ij = 0, or that overlap (overlaps) h ij < 0.
• Contacts play a crucial role in the stability of the system. Their number is fixed to be exactly equal to the number of degrees of freedom and its fluctuations are suppressed, as it happens at jamming [7]. We show that, as at jamming [18], the spatial fluctuations of the local connectivity of the contact network are hyperuniform. Conversely, the fluctuations of the number of overlaps follow central limit theorem and spatial fluctuations of the overlap network are only uniform.
• If we look at gap variables and we focus on strictly positive and negative gaps, we find that both distributions have a power law divergence for small argument (in absolute value). The power law exponents controlling the divergence appear to be the same -within numerical precision -for both distributions and very close to the one of positive gaps at jamming.
• Contacts can be associated with forces in the interval [0, f c ]. We measure the force empirical distribution and show that it displays two singular pseudogaps, close to zero and close to f c . The pseudogap exponents appear to depend on the packing fraction close to jamming. However, if we carefully separate the contribution of "bucklers", namely spheres that have d + 1 interacting spheres [13], from the bulk statistics, the pseudogaps are universal and characterized by the same exponents in the whole jammed phase, far from jamming. The values of the critical exponents appear to be the same -within numerical precision -as the one of small force distribution at jamming.
Isostaticity and critical behavior in the force and gap distributions have been shown to appear in the unsatisfiable phase of the spherical perceptron optimization problem with linear cost function, which is a mean field model for linear spheres [19]. The main result of the present work is that these properties appear to survive in a robust way when we go to finite dimension. This implies that jammed packings of linear spheres are characterized by diverging isostatic lengthscales and therefore are critical even far from jamming in the compressed phase.

Numerical simulations
The linear ramp v 1 (h) is a singular interparticle potential and therefore, both for the sake of the theoretical comprehension and to perform numerical simulations, it is very useful to smooth the singularity out and to define a differentiable -regularized potential between particles. From now on we set f c = 1 1 . We can define a regularized interparticle potential as which in the → 0 reduces to the linear ramp potential. The potential energy of a system of N spheres is therefore defined as H (x) = i<j v 1 (h ij ; ) and we want to study what happens in the → 0 limit. In Eq. (1), the nondifferentiable point in the origin of v 1 (h) has been regularized by an arch of parabola of curvature 1/ . Numerically, the regularization enables us to use gradient based routines. Theoretically, the model splits the degeneracy of forces at h = 0 in an interval of order , i.e. − /2 < h < /2; it also allows to properly define the Hessian controlling minima's stability and to argue in favor of isostaticity. We consider systems of N up to 4096 disks in dimension d = 2 and N up to 1024 spheres in dimension d = 3, inside a d−dimensional box of side-length L with periodic boundary conditions. The particles' radii R {i=1...N } are random uniformly distributed between the values 1−p and 1+p, with polydispersity p = 0.2. The side-length L of the box is set by the volume density is Euler gamma function. Starting from a random configuration of particles' positions, we minimize the energy of the system H (x) = i<j v 1 (h ij ; ) with the regularized interparticle interaction potential v 1 (h ij ; ) defined in eq. (1). The first minimization is run with = 10 −2 using FIRE minimization algorithm [20]. From the obtained configuration, we reduce by a factor 2 and repeat the minimization using an approximated conjugate-gradient method (we use the routine L-BFGS [21]). We repeat the procedure halving at each step and we stop at ∼ 10 −8 . Using more accurate minimizers, it is possible to access to lower values of . We check that for the final value of we have that min ij |h ij |, where min ij |h ij | is the smallest non-zero gap of the configurations we are looking for. With this procedure we meet the jamming transition at packing fraction φ 2d J 0.84 and φ 3d J 0.64. This procedure provides the set C of the C = |C| particles pairs µ = ij , with i < j, that are in contact (i.e. − /2 < h ij < /2) and the set O of the O = |O| particles pairs µ = ij that are overlapping (i.e. that have negative gaps h ij < − /2). Associated to the contact pairs, there are the scalar contact forces f µ = f ij that form a C−dimensional vector f = {f ij }, while overlapping spheres exchange forces of intensity 1, whose corresponding O−dimensional vector is simply 1 = [1, ..., 1]. The scalar contact forces f ij can be computed from the regularized potential of eq. 1 as f ij = |h ij − 2 |/ , implying 0 < f ij < 1. Introducing the matrices S and T , with dimensions C × N d and O × N d respectively, defined as S kα ij = (δ jk − δ ik )n α ij , with ij ∈ C, and we can compute the contact forces f ij also in an algebraic manner using the force-balance condition (2) Notice that the system self-organizes in a way that the forces solving the linear system (2) lie in the interval (0, 1). In Fig.1, we show an example of a configuration we obtain through the numerical procedure just described. In red and black we draw respectively the contact and overlap networks. In the following, we present data for d = 2 and N = 512 (unless otherwise specified). The data we got in d = 3 is qualitatively similar to the d = 2 case and therefore we report it in the appendix.

The jammed phase
In the jammed phase, for ϕ > ϕ J , particles overlap and therefore the numbers of contacts C and of overlaps O, the energy E and the pressure p are different from zero. In two and in three dimensions, in all the minima we found, once removed the rattlers 2 , C is equal to the isostatic value (N * − 1)d, where N * is the number of spheres which are not rattlers. 3 On the other hand, O, E and p are continuous at jamming, having defined the pressure p as p = V −1 i<j |r ij |f ij /d, where f ij is the force exchanged by the spheres i and j which can be f c in the case of spheres overlapping or it can be in the interval (0, 1) in the case of spheres in contact. For small momentum, the structure factor of the contact network decreases down to zero implying hyperuniformity in the fluctuations of connectivity. The exponent controlling the behavior of the structure factor appears to be close to ∼ 1.53 which is the same found at jamming [18]. On the contrary, the connectivity of the overlap network is not hyperuniform. These data are produced with system size N = 4096. This is shown in Fig.2 where we plot the isostaticity index c = C/[(N * − 1)d], for a 2d system. In the same figure, we plot the overlap fraction , and, in the inset, the energy per particle and the pressure. These quantities start from zero at jamming and monotonically increase as the packing fraction grows. However we do not attempt a characterization of their scaling behavior close to unjamming [10] living it for future work.
To characterize the networks of interaction, we study the fluctuations in the local contact number and overlap number. Following [18], we look at the local connectivity fluctuations of the networks of interactions by measuring the structure factors where c i represents the number of contacts or overlaps of particle i for the contact or overlap structure factors respectively, δc i = c i − c is its local fluctuations and the angular brackets represent average over different minima. We plot both structure factors in Fig.3   of the overlap number are normal and scale as V . This difference is a manifestation of the different role that contacts and overlaps play in the stability of the system. As the system is progressively compressed from the jamming point to higher densities, the networks self-organize keeping the number of contacts fixed while increasing the overlaps. As at regular jamming [18], fluctuations of contact numbers away from isostaticity are suppressed and controlled by an infinite lengthscale. We note that the structure factor of the contact network shrinks to zero with a power law that is close to what is observed at jamming [18]. We conclude by noting that while increasing the packing fraction, the fraction of overlaps displays a plateau around ϕ ∼ 1.5. We empirically observe that at the same point the overlap network seems to undergo to a kind of percolation transition whose nature and properties are left for future investigations.

Statistics of gaps
Having established that the system is isostatic, it is natural to turn the attention to the distribution of non-zero gap variables, which at jamming provides an important characterization of criticality. While at jamming all gaps are positive or zero, here we also have 'negative gaps', quantifying the overlaps between particles. Both the distributions of positive and negative gaps appear to be singular at small argument. In Fig. 5 we plot the cumulative distribution of both positive and negative gaps for several packing fractions beyond the jamming transition point. The small gap behavior is controlled in both cases by a power law. If we denote by g ± (h) the positive and negative gap distribution, we have that at small argument. The two exponents coincide within the errors, γ + ≈ γ − and their numerical value appear to be independent of density and equal to the one of positive gaps at jamming γ ± = γ J ≈ 0.41 . . . [15], as predicted by mean field theory [19].

Statistics of forces
Local minima contain an isostatic number of contacts to which we can associate contact forces and study their empirical distribution. Scalar contact forces are naturally defined in the interval [0, f c = 1] and we observe that, as soon as we enter in the jammed phase, their distribution develops two pseudogaps close to the edges f ∼ 0, 1 (see Fig. 4). We observe that the exponents controlling the two pseudogaps appear to depend smoothly on the density gaps CDF |h|/ |h| The positive and negative cumulative gap distribution. The y-axis of the negative gap is rescaled by an artificial factor 0.1 to improve the readability of the figure. We observe that both cumulative distribution appear to be described by the same power law exponent for small argument.
especially close to jamming and for f ∼ 0 + . However, following [13] we perform a statistical analysis in which we remove bucklers, namely spheres interacting with d + 1 spheres. The result of the analysis is plotted in Fig. 4 and we show that, independently from the packing fraction, the force distribution behaves as . is the critical exponent controlling small forces between hard spheres at jamming [15]. In Fig. 4, we also plot the cumulative distribution function of bucklers' forces close to f ∼ 0. Again we see a power law behavior controlled -within numerical precision-by the same power law exponents controlling bucklers at jamming of hard spheres [11,13]. Finally, we note that deep in the jammed phase, localized effects such as bucklers (but also rattlers) disappear (with the statistics we have access to) and we do not need to separate them from the statistics of forces to observe a critical power law with mean field exponent θ J .

Non-linear marginal stability of linear soft spheres
The exact solution of the perceptron optimization problem with a linear cost function [19] provides a comprehensive mean-field framework that predicts the main features discussed in this paper namely isostaticity, identity of exponents θ + = θ − , γ + = γ − , their numerical values and so on. This theory can be adapted straightforwardly to soft linear spheres in infinite dimension [22]. In the next section we complement the mean field theory with marginal stability arguments.

Local stability
The configurations of minima of the PEL at finite energy density contain overlapping particles. It is easy to understand that there should also be pure contacts. The forces corresponding to overlapping particles are constant in modulus (equal to f c = 1) and, without contacts, only very symmetric configurations of particles would be mechanically stable. In fact, more generically, a number of contacts less than d on a particle would require a highly symmetric configuration to be stabilized by only overlaps. The minimal number of contacts necessary to stabilize a single particle is therefore d, with a number of overlaps larger or equal to one (or with at least another contact). When we go from the jammed phase towards the jamming point, the number of overlaps vanishes and we recover that at jamming a number of contacts larger or equal to d + 1 is required to block a sphere. Particles with d + 1 interactions are prone to local excitations and are usually called bucklers.

The regularized Hessian and isostaticity
Local minima of the linear ramp potential are anharmonic due to the singularity in the pairwise interaction potential. However, one can consider the -regularized potential and look at the Hessian of local minima in this case. This is indeed well defined and reads Focusing on i = j, we have the first term, often called prestress, which vanishes at jamming while we call the second term the elastic term. Because of the regularization, we have the indicator function which is equal to one if A is true and zero otherwise. Notice that the Hessian receives contributions both from overlaps and contacts. Overlaps contribute just to the prestress. Contacts instead contribute both to the prestress, with a finite term (notice that (h − 2 )/ is actually the contact force), and to the elastic part with a term proportional to 1/ . This implies that for a variation of the position of the particles such that |δx i | , the energy stored in the elastic term is of order , and dominates the one stored in the prestress which is of order 2 . This is a crucial property, which is at the basis of isostaticity and the criticality of non-linear excitations in the compressed phase. Despite giving only a relatively small contribution, the contribution of the prestress term is important. In fact, as usual in repulsive sphere systems, this is a destabilizing term (it corresponds to a negative definite matrix) that, though small, would imply unstable directions if the elastic part is not full ranked. We conclude that the number of contacts should be at least isostatic so that the total matrix is positive definite and the minimum is stable.
The Hessian is therefore dominated by its isostatic random elastic part. Isostatic random matrices are gapless [4,[23][24][25][26][27] and characterized by an abundance of soft modes, their spectrum should behave as λ −1/2 at small argument, where λ represents the eigenvalues. We measure the spectrum of the elastic term matrix, namely the spectrum of lim →0 H ab ij . In Fig.6, we plot the corresponding density of states (DOS) with respect to the vibrational frequency ω = √ λ. Varying the density from φ = 0.85 to φ = 2.0, our numerical simulations are compatible with having a constant DOS for ω → 0. In the appendix we develop a mean field theory for such behavior supporting this numerical finding.

Non-linear excitations
Further information can be gained considering a non-linear stability analysis for the local minima of the PEL. The data we presented clearly shows that minima are isostatic configurations where the distributions of both positive and negative gaps display a power law behavior at small argument. At the same time the isostatic delta peak of marginally satisfied gaps is accompanied by a contact force distribution which has two pseudogaps close to forces equal to zero or one. The emergence of these power laws controls the non-linear excitations that dominate the dynamics of the system when perturbing it away from such local minima. One can understand the nature of those excitations generalizing the lines of reasoning employed in [11,14] for the jamming point. The simplest excitations are the ones in which isostaticity is off by one contact. There are here two possibilities, either separating two spheres in contact and opening a positive gap, or on the contrary pushing two spheres in contact to make them overlap and create a negative gap. The softest excitations are then the ones corresponding to either very week contacts in the former case, or to contacts carrying a force close to one in the latter case. When such contacts are removed, the system would become mechanically unstable unless a new contact forms in the system and again we have two possibilities, either a gap closes, or an overlap relaxes to become a contact. Assuming that both processes occur with finite probability, we have θ + = θ − and γ + = γ − . Following [11,14], one arrives at the scaling relation γ + = 1/(2 + θ + ) controlling the critical exponents, which is verified by both our numerics and the mean field theory of [19].

Discussion and Conclusion
The main conclusion of our analysis, both theoretical and numerical, is that linear soft spheres are critical even above the jamming transition, in the overcompressed phase. The criticality of local minima of the PEL is described by a set of power laws controlling the positive and negative gap distributions as well as contact forces. The critical exponents controlling such distributions appear to be numerically indistinguishable from the corresponding ones at jamming. Furthermore, the critical behavior is again directly controlled by isostaticity of local minima. This is an interesting result that opens the way to study jamming criticality in a different and complementary way. Indeed, typically, in order to look for the critical properties of the jamming transition, one needs to fine-tune the numerical simulations in order to be close to jamming. Linear soft spheres instead allow us to get to jamming-like critical configurations just by looking at local energy minima which can be obtained using standard numerical routines to minimize the energy. In this case no fine-tuning is needed. The rich physics that we observe in linear spheres is due to isostaticity which we robustly find with descent dynamics in local minima at finite N . Our work opens a series of perspectives: on one hand, it would be extremely interesting to characterize the rheology of strained linear spheres [28,29]. A possible way to look for that would be to perform similar experiments as in Ref. [30][31][32] and to analyze the statistical properties of avalanches. On the other hand, it would be interesting to investigate other concave penalty exponents α < 1, or more complex potentials, to see if different non-linear criticality may arise. Moreover, by switching on temperature, one may investigate if marginal stability emerges at a critical point, the Gardner transition [33][34][35]. Finally, further work is required to understand the behavior of bulk quantities such as energy and pressure close to unjamming. Likely, this cannot be obtained from the scaling valid for α > 1 [4,10,36,37], and it may be important to study the leading corrections to this scaling for α close to one. While our data hints at such phenomenology, further investigations are needed. A possible way to investigate this point would be to progressively compress a configuration from jamming. The dynamics should be dominated by contacts carrying forces close to one becoming overlaps while small gaps becoming contacts with a net flux of gaps from the positive to the negative side of the distribution. How to describe such dynamics is left for future work. Fig. 7-left we plot the dependence of the fraction of contacts with respect to degrees of freedom (computed removing rattlers). In the same plot we also plot the fraction of overlaps. While we see that the system is isostatic at all packing fraction above jamming, the density of overlaps increases monotonically. In the inset we report the behavior of the energy and pressure as a function of φ. In the right panel we plot also the structure factor of the local connectivity of the overlap and contact network. We show that while the overlap network obeys central limit theorem, the fluctuations of the number of contacts are suppressed and the structure factor goes to zero for small momenta. Therefore both in d = 2 and d = 3 isostaticity is reached when minimizing the energy of the system. In order to see how this happens numerically, in Fig. 8 we plot the isostaticity index as a function of the regularization parameter (we plot data for d = 2 for simplicity). We clearly see that as soon as the linear potential limit is reached, the system self organizes to sit on an isostatic minimum. Finally in Fig.9 we plot the contact force and gap distribution. Contact forces display two pseudogaps close to the edges f = 0, 1 while small gaps are controlled by a power law distribution. The critical exponents controlling these distributions appear to be the same (within numerical precision) with the ones of two dimensional linear spheres.
Therefore the main conclusion of this analysis is that, as for the jamming transition, the properties of soft spheres interacting with linear potential do not depend on the dimensionality of the system (apart from local bucklers/rattlers effects). gaps CDF |h|/ |h| In the main text we have argued that the elastic part of the contact network matrix has a density of states D(ω) which goes to a positive constant for ω → 0 analogously to what happens at jamming. In this section we construct the mean field theory for such behavior. We consider the spherical perceptron optimization problem with linear cost function studied recently in [19], see also [38][39][40]. This model is in the same universality class of linear soft spheres. It is an optimization problem where a set of N variables x i arranged in a vector x = {x 1 , . . . , x N } lying on the sphere |x| 2 = N are sought to minimize the cost function where the gaps h µ are defined as with ξ µ = {ξ µ 1 , . . . ξ µ N } a set of N dimensional vectors with components extracted from a Normal distribution and σ a constant control parameter. Local minima of the PEL are isostatic meaning that there is an isostatic number of gaps h µ = 0 and characterized by critical power laws in the gap and forces distribution. The Hessian of the non-analytic minima can be defined by smoothing out the singularity of the linear potential close to h µ = 0 as we have done with linear soft spheres. Calling the smoothing parameter, the Hessian becomes being ζ a Lagrange multiplier needed to enforce the spherical constraint which plays here the same role of the prestress in spheres. In the glassy phase, ζ < 0 and therefore for → 0 one needs to have isostatic minima. Since the system is isostatic, assuming that the patterns are random and the only ingredient that matters for the statistics of the Hessian is the number of contacts [26], we get that the spectrum of the Hessian is a given by a Marcenko-Pastur distribution for the eigenfrequencies ω given by This result holds in the whole glassy phase regardless of the energy, being the glassy phase always isostatic. In Fig. 10 we plot the density of states as extracted from numerical simulation for α = 5 and at two different values of the energy and we compare it with the theoretical prediction.