Jamming and replica symmetry breaking of weakly disordered crystals

We discuss the physics of crystals with small polydispersity near the jamming transition point. For this purpose, we introduce an effective single-particle model taking into account the nearest neighbor structure of crystals. The model can be solved analytically by using the replica method in the limit of large dimensions. In the absence of polydispersity, the replica symmetric solution is stable until the jamming transition point, which leads to the standard scaling of perfect crystals. On the contrary, for ﬁnite polydispersity, the model undergoes the full replica symmetry breaking (RSB) transition before the jamming transition point. In the RSB phase, the model exhibits the same scaling as amorphous solids near the jamming transition point. These results are fully consistent with the recent numerical simulations of crystals with polydispersity. The simplicity of the model also allows us to derive the scaling behavior of the vibrational density of states that can be tested in future experiments and numerical simulations.


I. INTRODUCTION
Physics of crystal and amorphous solids are qualitatively different.For instance, low frequency eigenmodes of crystals are phonon, and thus the vibrational density of states D(ω) follows the Debye law D(ω) ∼ ω d−1 , where d denotes the spatial dimensions [1].On the contrary, amorphous solids have excess nonphonon excitations.As a consequence, the density of states normalized by the Debye's prediction D(ω)/ω d −1  shows a peak at a certain frequency ω = ω BP [2][3][4][5][6].This phenomenon is known as the boson peak and thought to be one of the universal properties of amorphous solids [7].
Crystal and amorphous solids also show distinct elastic properties near the (un) jamming transition point at which constituent particles lose contact, and simultaneously the pressure vanishes [8].Here we focus on the jamming of spherical and frictionless particles interacting with finite and repulsive potentials.The scaling of these models is now well understood due to extensive numerical simulations [8,9] and theories [10][11][12][13].The shear modulus G of crystals does not show the strong pressure p dependence and remains a constant at the jamming transition point [14,15].On the contrary, G of amorphous solids shows the power law behavior G ∼ p 1/2 and vanishes at the jamming transition point [8,9].The behavior of G is directly related to the contact number per particle Z as G ∝ δZ ≡ Z − Z iso [11].Here Z iso denotes the contact number when a system is isostatic, i.e., the number of constraints is the same as the number of degrees of freedom [16,17].At the jamming transition point, δZ > 0 for perfect crystals, leading Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
Crystal and amorphous are two extreme states of solids: the former is a state free from disorder while the latter is a state of maximum disorder.From both theoretical and practical points of views, it is important to understand how the physical properties shift from that of crystal to amorphous on the increase of the strength of disorder.Previous numerical simulations show that small disorder only play a moderate role far from the jamming transition point (p ∼ 1).For instance, numerical studies of crystals with polydispersity show that the amplitude of the boson peak D(ω BP )/ω d−1 BP only continuously increases on the increase of the polydispersity η, if η is small enough [18,19].Near the jamming transition point p 1, on the contrary, even small disorder dramatically change the physical properties of crystals.More and more nonphonon modes appear as p decreases, eventually leading to the divergence of D(ω BP )/ω d−1 BP in the jamming limit p → 0, in sharp contrast to perfect crystals where D(ω) does not show the strong p dependence.Furthermore, for crystals with small defects or polydispersity, G and δZ exhibit the same power laws of amorphous solids sufficiently near the jamming transition point [14,15].In particular, G and δZ vanish at the jamming transition point if there is even infinitesimally small polydispersity [15,20], while for perfect crystals, G and δZ remain finite.
Our aim here is to construct a solvable mean-field model being able to describe the above striking effects of disorder on crystals near the jamming transition point.We consider a model in the limit of large dimensions, which is a popular mean-field limit in theoretical physics [21,22].In this limit, only the first virial corrections give a relevant contribution [22,23].For a short-range potential such as hard spheres, this implies that the information of nearest neighbor structures is enough to describe the physics.Motivated by this consideration, we introduce an effective single-particle model that only takes into account the interactions between a particle of interest and nearest-neighbor particles.For zero polydispersity η = 0, our model correctly reproduces the scaling of perfect crystals.For finite η, on the contrary, our model predicts that the existence of the replica symmetry breaking (RSB) transition [24] at finite pressure p = p RSB .For p p RSB , the model exhibits the same scaling as amorphous solids.Thereby, our model can reproduce the sharp cross-over from the scaling of crystal to amorphous observed in previous numerical simulations of weakly disordered crystals.

II. MODEL
We consider a tracer particle surrounded by M frozen nearest-neighbor (NN) particles in the limit of the large spatial dimensions.In d spatial dimensions, the tracer particle has d degrees of freedom, implying that the model becomes isostatic when the contact number is Z = Z iso = d. 1 The tracer and NN particles interact with the one-sided harmonic potential: where v(x) = x2 θ (−x)/2, and the prefactor denote the positions of the tracer and μth NN particle, respectively.σ μ denotes the interaction range between the tracer and μth NN particle.We assume that σ μ can be written as where σ controls the mean size of particles, and b μ denotes the polydispersity.The pre-factor of b μ , 1/d, is necessary to keep the relative interaction volume remains finite in the large dimensional limit: lim d→∞ (σ μ /σ ν ) d = e b μ −b ν .b μ follows the normal distribution of zero mean and variance η 2 : The NN particles are homogeneously distributed on the surface of the hypersphere of the radius √ d, namely, the distribution function of y μ is given by To get the physical intuition about the model, we first explain the behavior at the jamming transition point in d = 2 in the absence of the polydispersity η = 0, comparing it with the hexagonal packing.For the hexagonal packing, the NN particles are arranged periodically on the equidistant line from the tracer particle [Fig.1(a)].On the contrary, for our model at the jamming transition point, the NN particles are randomly distributed on the equidistant line [Fig.1(b)].The tracer is in contact with all NN particles, as in the case of hexagonal packing, leading to a hyperstatic configuration when the number of the NN particles M is larger than Z iso .
The same story holds in general d as long as η = 0: the jamming occurs when σ ≡ σ 0 J = √ d, at which X = 0 and h μ = 0 for all μ, meaning that the tracer particle is in contact with all NN particles, leading to a hyperstatic configuration when M > Z iso .On the contrary, for η > 0, the jamming configuration is nontrivial, which we shall discuss in this manuscript.

III. MARGINAL STABILITY
The previous works for the mean-field models of the jamming transition unveiled that the systems undergo replica symmetry breaking (RSB) before reaching the jamming transition point.In the RSB phase, the systems are marginally stable [24,25].Here we show that the contact number in the RSB phase can be calculated by using the marginal stability.
At zero temperature, the stability of the system can be discussed by observing the Hessian of the interaction potential: where the 1/d prefactor is necessary to make H i j = O(1), and we only keep the relevant terms in the limit of d → ∞.In this limit, y μ i can be identified with the i.i.d.Gaussian random variable of zero mean and unit variance, see Appendix B. Thus, H i j can be considered as a Wishart matrix [26] with an additional diagonal term.The eigenvalue distribution of H i j follows the Marchenko-Pastur distribution [27]: where we have defined the contact number per degree of freedom z = Z/d and pressure p as In the RSB phase, the marginal stability requires λ − = 0 [24,25].This condition with Eq. ( 6) determines z as a function of p: This result implies that (i) the model is isostatic z = z iso = 1 at the jamming transition point p = 0 2 and (ii) z exhibits the square root scaling δz = z − z iso ∝ p 1/2 for p 1.Those properties are the same as amorphous solids consisting of soft harmonic particles [9] and a mean-field model of the jamming transition [27].
Here we used a rather heuristic argument to calculate z in the RSB phase.But the same result can be derived by directly solving the RSB equation, see Appendix E.

IV. REPLICA METHOD
The RSB is a consequence of the complex structure of the free energy landscape of amorphous solids near the jamming transition point [13].On the contrary, perfect crystals or nearly perfect crystals have a unique minimum, in other words, the systems are in the replica symmetric (RS) phase.Since the RS phase is not marginally stable, we can not use Eq. ( 7).Instead, we calculate z by using the replica method [24].The calculation is very similar to that of the perceptron, which was previously investigated as a mean-field model of the jamming transition [27,28].Therefore below we just briefly sketch how to calculate z as a function of p.The details of the calculations are provided in Appendix C.
To calculate z and p, it is convenient to introduce the gap distribution: where F denotes the free energy, and • denotes the average for both quenched disorder and thermal fluctuation.Using g(h), z and p are calculated as We calculate F by using the replica method: where Z denotes the partition function: β represents the inverse temperature.In this work, we investigate the model only in the zero temperature limit β → ∞.
Investigations at finite β are left for future work.The square brackets in Eq. ( 11) denote the average for the quenched randomnesses: 2 The condition of the isostaticity is now z = z iso ≡ Z iso /d = 1.where P(b μ ) and P(y μ ) are given by Eq. ( 3) and ( 4), respectively.By using the saddle point method, one can represent F as a function of the correlation among the replicas (see Appendix B for details): where X a and X b denote the positions of the ath and bth replicas, respectively.As the free energy has a single minimum in the RS phase, there is no reason to distinguish a specific pair of replicas ab, implying that Q ab is written as We can calculate q and q 0 by solving the saddle point equations: ∂ q 0 F = 0 and ∂ q F = 0 (see Appendix C for details).Substituting back the results to F , we obtain the free energy at the saddle point, which allows us to calculate g(h), z, and p. Below, we will show the results only for M = 10d, but we confirmed that the qualitatively same results are obtained for different values of M.
To see the stability of the RS ansatz, we calculate the minimal eigenvalue λ − by substituting the RS results for z and p into Eq.( 6).The RS-RSB transition point is determined by the condition λ − = 0.In Fig. 2, we show the RS-RSB phase diagram in the η − p plane.It is noteworthy that the RSB always occurs at a finite pressure p = p RSB ∼ η 2 before reaching the jamming transition point p = 0 whenever η > 0.

V. SCALING OF CONTACT NUMBER
Following the above procedures, we calculate z for several η.We summarize our results in Fig. 3.There are three different scaling regions.For p η, z takes a constant value z ≈ M/d, meaning that the tracer particle contact with most NN particles, see the black line.For η 2 p η, the contact number decreases as δz ∼ p, see the blue dotted line.At p = p RSB ∼ η 2 , the RS solution becomes unstable, and for p p RSB , one should use the RSB result Eq. (8).For p η 2 , Eq. ( 8) predicts δz ∼ p 1/2 , see the red dashed line.
For p p RSB , the results for different η collapse on a single curve if one plots δz as a function of η −1 p, see Fig. 3(b).This scaling is consistent with a previous numerical simulation [15] and perturbation theory [29].Remarkably, the above scaling implies that the two limits η → 0 and p → 0 are not commutative: if one takes the limit η → 0 first and then takes the limit p → 0, one gets δz > 0, contrary, if one takes the limits in reverse order, one gets δz = 0.

VI. DENSITY OF STATES
An important quantity to characterize the physics of solids is the vibrational density of states D(ω), which is a distribution of the eigen-frequency ω = √ λ.By using Eq. ( 6), D(ω) is calculated as D(ω) = 2ωρ(ω 2 ).Near the jamming transition point for small ω, D(ω) asymptotically behaves as where ω 0 = √ λ − .In the RS phase, ω 0 > 0 and D(ω) has a finite gap. 3 ω 0 decreases on the decreasing of p and eventually vanishes at p = p RSB .In the RSB phase p p RSB , ω 0 = 0 and D(ω) is gapless.For ω 1, the density of states exhibits the quadratic scaling D(ω) ∼ ω 2 .This is the same result as previous mean-field theories of amorphous solids [25,30].In the jamming limit p → 0 for η > 0, D(ω) always exhibits the plateau for small ω, which is fully consistent with previous numerical simulations of weakly disordered crystals near the jamming transition point [14,31].Now we want to calculate the boson peak.For comparison with numerical simulations, we consider the height of D(ω)/ω m at its peak ω = ω BP , where m = 1 and m = 2 correspond to the Debye predictions in two and three spatial dimensions, respectively.Using the scaling of δz and ( 16), one can deduce the asymptotic behavior for m 24 as a function of p: Equation (17) suggests that the boson peak intensity diverges in the jamming limit p → 0. This scaling is the same of that of amorphous solids near the jamming transition point observed by a numerical simulation of three dimensional harmonic spheres [32].Repeating the similar calculation, one can derive the scaling of the boson peak intensity as a function of η: Equation ( 18) suggests that, on the increase of the polydispersity η, the boson peak begins to increase at η ∼ p.This is consistent with a previous numerical simulation of crystals with small polydispersity [15].In Fig. 4, we summarize the scaling of the boson peak intensity predicted by the above equations.It is interesting to test the full scaling behavior by experiments and numerical simulations.

VII. SUMMARY AND DISCUSSION
In this work, we have introduced a mean-field model to describe the jamming transition of crystals with small polydispersity.We solved the model by using the replica method and determined the full scaling behaviors of the contact number and density of states above the jamming transition point.The results are well agreed with previous numerical simulations.
Another important quantity to characterize the jamming transition is the gap distribution [12].In Refs.[31,33], it is shown that the gap distribution of the disordered crystal has a different critical exponent from both perfect crystals and amorphous solids.It is an interesting future work to see if our model can explain this intriguing behavior of the gap distribution.

APPENDIX A: INTERACTION POTENTIAL IN THE LARGE DIMENSIONAL LIMIT
Our model consists of a tracer particle and M nearestneighbor (NN) particles.The tracer is located in a ddimensional hypersphere, and the M NN particles are fixed on the surface of the hypersphere.The interaction potential is given by where v(h) = h 2 θ (−h)/2 and Here the prefactor √ d is necessary to keep V (X ) = O(d 0 ) in the d → ∞ limit, as we will see below.X = {X 1 , . . ., X d } and y μ = {y μ 1 , . . ., y μ d } denote the positions of the tracer and μth NN, respectively.The distribution of y μ is σ μ denotes the interaction range between the tracer and μth obstacles.We consider that σ μ has the following form: where σ controls the mean size of the particles, and b μ denotes the polydispersity.b μ follows the normal distribution of zero mean and variance η 2 : For η = 0, the jamming transition occurs at σ = σ 0 J ≡ √ d at which X = 0 and hμ = 0 for all μ.We expand σ around σ 0 where the pre-factor of a, 1/d, is necessary to keep the relative interaction volume (σ /σ 0 J ) d finite in the limit of d → ∞.Substituting Eqs.(A4) and (A6) into the gap function hμ and expanding by 1/d, we get where we have used y μ • y μ = d.We require that the firstorder terms have the same magnitude.This is possible if the following conditions are satisfied ).We introduce a new variable of order one: Equations ( A9) and (A10) lead to d i=1 x i y μ i = O( √ d ), which is a natural result because y μ i is a random variable of zero mean and unit variance.Up to the first order, we get the following result: In summary, in the limit of d → ∞, the interaction potential is The gap function h μ has a similar form to the perceptron, except the additional terms x • x/2d and b μ [27,28].Therefore we can apply the same technique of that of the perceptron to investigate the model.

APPENDIX B: FREE ENERGY
Although we only investigate the model at zero temperature T = 0, we fist consider the free-energy at finite T to apply the technique of the statistical mechanics and then take the limit of T → 0. The free energy can be written as where and Here we introduced the inverse temperature β = 1/T .First, we will show that the average for y μ i can be replaced by the average for a normal distribution of zero mean and unit variance.The mean value of f (y μ ) = ln Z is represented as In the limit d → ∞, we can evaluate the integral of λ by using the saddle point method.The saddle point condition is where we used y μ • y μ = O(d ) and ln f = ln ln Z = O(ln d ) d. Applying the saddle point method for the integral of λ in Eq. (B4), we get meaning that the distribution function of y μ i converges to a normal distribution of mean zero and unit variance in the limit of d → ∞.This can greatly simplify the calculation as we will see below.Now, we calculate the free energy Eq.(B1) by using the replica method: Using the Fourier transformation of the delta function δ(x) = (2π ) −1 d re irx , the partition functions can be written as where a = 1, . . ., n and μ = 1, . . ., M. Since y μ and b μ follow the normal distribution, one can show that where we have introduced a new variable The Jacobian of the change of the variables is where and Q * ab detenos the saddle-point value satisfying Finally, the free-energy Eq. (B7) is calculated as APPENDIX C: CALCULATION WITH THE REPLICA SYMMETRIC ANSATZ

Free energy
Here we investigate the model by assuming the replica symmetric (RS) ansatz: Then, the free-energy is where we used the abbreviations: For low T , we can perform the harmonic expansion: Substituting it into the free-energy, Eq. (C2), we get in the T → 0 limit χ and q 0 are determined by the saddle point equations: (C6) The equations can be solved numerically, which allows us to calculate q 0 and χ for given a and η 2 .

Gap distribution function
Our goal is to calculate the contact number z as a function of the pressure p.For this purpose, it is convenient to introduce the gap distribution function: where • denotes the average for both thermal fluctuation and quenched disorder.g(h) can be calculated as In the limit of T → 0, the saddle point method leads to The contact number per degree of freedom z and pressure p are calculated from g(h) as

Numerics
We calculate z as a function of p with the following steps.(1) Calculate q 0 and χ as functions of a and η by solving Eqs.(C6).( 2) Calculate z and p by substituting the above results into Eq.(C10) and (C11).( 3) Plot z as a function of p.
We found that the above algorithm does not converge for M ∼ d.This may imply that if the number of the NN particles is not large enough, the tracer particle can escape, and the harmonic expansion, Eq. (C4), breaks down.To avoid this problem, in the main text, we show the results for M = 10d d.We checked that the qualitatively same results are obtained for different values of M as long as M d.

Scaling
We derive the scaling behavior near the jamming transition point for η 1 from the asymptotics of the RS equations.First, we discuss the scaling at the jamming transition point.At the transition point, the harmonic expansion breaks down, meaning that χ → ∞.Therefore Eqs.(C6) reduce to By solving the above equations, one can calculate a and q 0 at jamming for given η.From an asymptotic analysis for η 1, we can show that The results is consistent with a naive dimensional analysis: both a and η have the dimension of length, leading to a ∝ η, and q 0 has the dimension of the squared of length, leading to q 0 ∝ a 2 ∝ η 2 .
To get the scaling above jamming, we rewrite the saddle point equations Eqs.(C6) as Using Eq. (C14), we get where we used Eq. ( C13) and e RS ∼ h 2 ∼ h 2 ∼ p 2 .Substituting it into Eq.(C15), we get where c 1 and c 2 denote constants.On the contrary, far from the jamming transition point, we get z ∼ α as the tracer contact with most NN particles.Summarizing the results, the RS ansatz predicts the following scaling: (C18) Note however that the above result for p η 2 is incorrect because the RS solution becomes unstable.As discussed below, the correct scaling δz ∼ p 1/2 is obtained by using the RSB equations.As discussed in the main text, if one plots z − 1 as a function of η −1 p, the results for p η 2 collapse on a single master curve.

APPENDIX D: EIGENVALUE DISTRIBUTION
We consider the Hessian matrix: where we have dropped the subleading terms in the large dimensional limit.From the central limit theorem, the diagonal term converges to the pressure: In the previous sections, we have shown that y μ i follows the normal distribution of zero mean and unit variance.Therefore M i j is a Wishart matrix shifted by p [26].The eigenvalue distribution follows the MarchenkoPastur (MP) law [25]: In particular, we are interested in the minimal eigenvalue λ − = ( √ z − 1) 2 − p, which vanishes when the RS Ansatz becomes unstable, and the RSB phase appears.transition point.By using the scaling in the RS phase, Eq. (C18), we can see that λ − for η 2 p η behaves as where c 1 denotes a positive constant, meaning that the RSB occurs at This is consistent with the numerical solution of the RS equations presented in the main text.

APPENDIX E: FULL RSB ANALYSIS
Here we calculate the contact number z as a function of p in the RSB phase by directly analyzing the full RSB free energy.

FIG. 1 .
FIG. 1. Jamming configurations for d = 2 and η = 0.The gray disk denotes the tracer particle, and the black disks denote nearest neighbors.The blue solid line denotes the equidistant line from the tracer particle, and the red dashed lines denote contacts.(a) Configuration of the hexagonal close-packing.The nearest neighbors are arranged periodically on the equidistant line.(b) Configuration of our model.The nearest neighbors are placed randomly on the equidistant line.

FIG. 3 .
FIG. 3. Scaling of the excess contact number δz as a function of the pressure p for M = 10d.(a) Filled circles denote the exact results.denotes the RSB transition point.The black solid line denotes z = M/d, the blue dotted line denotes δz ∼ p, and the red dashed line denotes δz ∼ p 1/2 .(b) The same data with the rescaled pressure.

)FIG. 4 .
FIG. 4. Scaling of the boson peak.(a) and (b) show the p and η dependence of the boson peak intensities, respectively.