Phase transition in the binary mixture of jammed particles with large size dispersity

It has been well established that particulate systems show the jamming transition and critical scaling behaviors associated with it. However, our knowledge is limited to (nearly) monodisperse systems. Recently, a binary mixture of jammed particles with large size dispersity was studied, and it was suggested that two distinct jammed phases appeared. Here, we conduct a thorough numerical study on this system with a special focus on the statistics of and finite-size effects on the fraction of small particles that participate in the rigid network. We present strong evidence that two distinct jammed phases appear depending on the pressure and composition of two species, which are separated by the first-order phase transition. In one of two phases, only large particles are jammed, whereas both small and large particles are jammed in the other phase. We also describe the phase diagram in the pressure-composition plane, where the first-order phase transition line terminates at a critical point. In addition, we investigate the mechanical properties in terms of the elastic moduli over the phase diagram and find that discontinuous changes in elastic moduli emerge across the phase transition. Remarkably, despite the discontinuities, the elastic moduli in each jammed phase exhibit identical scaling laws to those in the monodisperse systems.

It has been well established that particulate systems show the jamming transition and critical scaling behaviors associated with it. However, our knowledge is limited to (nearly) monodisperse systems. Recently, a binary mixture of jammed particles with large size dispersity was studied, and it was suggested that two distinct jammed phases appeared. Here, we conduct a thorough numerical study on this system with a special focus on the statistics of and finite-size effects on the fraction of small particles that participate in the rigid network. We present strong evidence that two distinct jammed phases appear depending on the pressure and composition of two species, which are separated by the first-order phase transition. In one of two phases, only large particles are jammed, whereas both small and large particles are jammed in the other phase. We also describe the phase diagram in the pressure-composition plane, where the first-order phase transition line terminates at a critical point. In addition, we investigate the mechanical properties in terms of the elastic moduli over the phase diagram and find that discontinuous changes in elastic moduli emerge across the phase transition. Remarkably, despite the discontinuities, the elastic moduli in each jammed phase exhibit identical scaling laws to those in the monodisperse systems.

I. INTRODUCTION
Jammed particulate systems are ubiquitous in our lives. Emulsions, colloidal suspensions, and granular materials are examples of jammed systems, where constituent particles are randomly jammed [1]. The random structure of jammed systems is one source of difficulties in understanding these systems.
One of the simplest models for jammed particulate systems is the assemblies of athermal particles that interact via short-range repulsive potentials. When we compress these particles from the low density, the system gains rigidity at the density called the jamming point. This phenomenon is known as the jamming transition established by many previous works, e.g., Ref. [2]. The geometrical, mechanical, and vibrational properties of the system are known to follow the critical power-law near the transition where the distance from the jamming points plays the role of a control parameter [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. Recently, jammed systems composed of dimer-shaped particles were shown to exhibit the same critical laws as a sphere packings [17][18][19]. In the last two decades, numerical and theoretical studies of these systems have developed to a large extent, and we have a good level of understanding of the critical behaviors near the jamming transition in nearly monodisperse jammed particles.
However, the jammed particles in reality are not monodisperse and are composed of particles of multiple sizes [20][21][22]. In particular, when the sizes of larger and smaller particles are quite different, these systems can exhibit phenomena that are never observed in the monodisperse systems [23,24]. For example, mixtures * hara-yusuke729@g.ecc.u-tokyo.ac.jp of large colloidal particles and small polymers are known to show rich phenomena, including the emergence of two distinctive colloidal glass phases and a peculiar nonlinear mechanical response [25][26][27][28].
Despite its theoretical and practical importance, the impacts of the polydispersity with large size ratio on the characteristic features of jammed packings are poorly understood. The simplest system to study the effects of polydispersity can be a binary mixture of particles with large size dispersity. The packings of binary mixtures have been studied to improve the packing efficiency. It is now well known that an injection of fine particles to coarse large structures realizes the random close packing with a higher packing fraction than monodispersity [29][30][31][32][33][34][35].
In addition, the nature of the jamming transition in binary systems was studied. Xu et al. showed that binary mixtures of equal numbers of small and large particles did not change the critical behavior near the jamming transition even if their sizes were significantly different [36]. Kumar et al. studied the effects on the bulk moduli by injecting a few percent of fine particles into a volume of large particles [37]. More recently, it was suggested that this system exhibited two different types of jammed phases, and there was a transition between them [34,38]. One phase is characterized by the jamming of only large components, and the other is characterized by the jamming of both components. Two different types of jammed states are also observed for two dimensional discs [33].
In the present work, we further investigate the structural and mechanical properties of the packing of binary mixtures of particles with large size dispersity. In particular, we thoroughly study the transition between two jammed phases, which was very recently reported in Refs. [34,38]. We focus on the fraction of small particles that contribute to the rigidity of the system as the order parameter and study their statistics and finite-size effect. We show strong evidence that there is a phase transition between two phases and that the transition is first-order in nature. We also reveal the phase diagram of this system and find that the first-order transition line terminates at a critical point. In addition, we show that the mechanical properties exhibit a discontinuous change at the first-order transition.
This paper is organized as follows. In Sec. II, we introduce the detailed description of the system in interest. In Sec. III, we present our simulation results. In Sec. III A, we focus on the analysis of the structural characterization at the fixed pressure and establish the first-order transition. Then, in Sec. III B, we extend the analysis to a broad range of pressure and reveal the phase diagram. Finally, in Sec. III C, the mechanical properties such as the bulk and shear moduli are studied. We provide the summary and conclusion in Sec. IV.

A. System description
We study the jammed packings composed of the binary mixtures of large and small particles in three spatial dimensions. The diameter of the large particle is 6 times that of the small particles. We denote the volumes of large and small particles as v l and v s , respectively, where v l = 6 3 v s . The numbers of large and small particles are denoted by N L and N S , respectively. These particles interact via a finite-range, purely repulsive potential [2]; two particles interact with each other only when they are in contact. The interaction potential is the following harmonic potential: where σ ij is the sum of the radii of particles i and j, r ij is the distance between the centers of particles i and j, and Θ(x) is Heaviside step function; Θ(x) = 1 for x > 0, otherwise Θ(x) = 0. We set the unit of the length and the energy as the diameter of the small particle and , respectively. The state point of this system can be characterized by two relevant parameters. The first one is the volume ratio of small particles: In this work, we focus on the range of 0.05 ≤ X S ≤ 0.25, where the number of small particles is approximately 10 to 70 times that of large particles. The second is the pressure: where ij is the summation over all contacting pairs of particles. The packing fraction is sometimes used as the control parameter, but we use the pressure because it enables us to more properly study the critical behaviors. Packings with desired (P, X S ) are obtained by iterations of compression and decompression of the system, as described below.
There are three possible quantities to define the system size: N S + N L , N S , and N L . Although previous studies [34,38] of binary mixtures use the total number of particles N S + N L to express the system size, we can also use the number of small particles N S or large particles N L . In this work, we use N S to precisely study the size effect on the susceptibility (see Sec. III A), while we use N L to explore the broad range of X S (see Sec. III B).

B. System preparation
We generate the jammed packing of desired pressures using the iterative compression and decompression of the system and the FIRE algorithm to minimize the energy [39]. The condition to terminate the relaxation al- Greek index α is the spatial index, and Roman index i is the index of the particles.
We prepare the configuration below the jamming point by relaxing the random configuration with a low packing fraction to the mechanical equilibrium. Then, we compress the system with fixed increment ∆φ ini = 0.01, and the systems are relaxed to mechanical equilibrium in each step. We continue to increase the packing fraction until the pressure of the equilibrated configuration exceeds the target pressures. The procedures are terminated if the pressure of the configuration is consistent with the target pressures within the acceptable error; in this case, δP = |P − P tag | ≤ 10 −2 P tag . When this condition is not satisfied, we (1) decompress with the rate ∆φ and (2) compress with the new rate ∆φ new = 0.5∆φ. The iterative use of the procedures brings the system to the desired pressure.
We generate at least 100 packings at each state point. We denote the average of physical quantities x in the ensemble of the packings generated at a state point by x sample . We also analyze the statistics of quantities x in the ensemble and calculate the probability distribution denoted as P (x). Note that the ordering of large particles occurs in some packings. We detect these packings by calculating the local bond order parameter of the structure of large particles [40]. These packings are rare, and we exclude them from our statistical ensemble.

A. Transition between two jammed phases
In this section, we study the structural properties of the packings with changing X S at the fixed pressure P = 10 −3 . We will show that the system exhibits two distinct jammed phases, which are separated by a first-order phase transition. One phase is named the L phase, where only large particles are jammed, and the other is the LS phase, where both large and small particles are jammed. To characterize these phases, we focus on the fraction of small particles that participate in the connected network of particles: Here, N r S is the number of small particles that are rattlers, which are defined as the particles whose contact number is less than d + 1 = 4 (where d = 3 is spatial dimension). We prepare at least 200 packings at each state point (P = 10 −3 , X S ) and calculate R s for each packing.
First, we focus on R S sample , which is the mean value of R S in the ensemble of at least 200 packings at each state point. Fig. 1 shows R S sample versus X s for various system sizes N S . R S sample monotonically increases with X s ; namely, a larger fraction of small particles corresponds to more small particles participating in the connected network. Notably, the change in R S sample becomes steeper when the system size increases. This implies there is a phase transition between the states with R S ≈ 0 and R S ≈ 1 in the thermodynamic limit, and R S is an order parameter in this phase transition.
To clarify whether this is a genuine phase transition, we analyze the statistics of R S . First, Fig. 2 shows the probability distribution of the order parameter R S over 400 different packings for several state points. When X S = 0.196, the distribution has a single peak at R s ≈ 0; most small particles are rattlers. When X S = 0.209, the distribution also has a single peak, but it is located at R s ≈ 1; most small particles participate in the connected network of particles. At X S = 0.202, the distribution shows the double peaks at R s ≈ 0 and 1, respectively; most small particles are rattlers in some packings, while most small particles participate in the connected network of particles in the other packings. The appearance of the bimodal distribution strongly suggests that this is the first-order phase transition between two phases. Next, we consider the system size dependence of this behavior. To facilitate the analysis, we introduce "susceptibility" χ, which is defined as Fig. 3(a) shows χ versus X S . Clearly, χ exhibits a peak, and the peak height increases with the system size. When the transition is first-order in nature, the peak height is expected to linearly depend on the system size. To study this point, Fig. 3(b) shows the system size dependence of the peak height of the susceptibility, which is denoted by χ peak . The data are consistent with the linear dependence on system size N S , which confirms that the transition between two phases is first-order. Finally, we evaluate the reduced cumulant of the order parameter defined by which is known as the binder parameter [41]. Fig. 4 shows U N S versus X S at various system sizes. Regardless of the system size, the binder parameter converges to 2 3 , which indicates the appearance of the bimodal distribution and a common feature of the first-order transition. This plot also shows the negative dip of the binder parameter, which is another common feature of first-order transitions. In summary, we provided strong evidences for the firstorder transition in the jammed phase, which is the transition between the jammed phase with only large components jamming (L phase) and that with both components jamming (LS phase).

B. Phase diagram
In the previous section, we established the presence of the first-order transition at P = 10 −3 , where fraction R S of small particles participating in the connected network is the order parameter. In this section, we extend the analysis to the broad range of pressures and determine the phase diagram of the packing of a binary mixture of particles. As a consequence, we will show that the first-order transition line terminates at a finite pressure. Fig. 5 plots R S sample versus X S at various pressure. At low pressure, R S exhibits discontinuous jumps around X S ≈ 0.2, which suggests the first-order transition as discussed in the previous section. However, at high pressure, e.g., P = 2 × 10 −2 , the increase becomes much milder. At this pressure, R S appears to smoothly change along X S , which implies the disappearance of the first-order transition.
To confirm this point in more detail, we analyze the statistics of R S at high pressure as performed at low pressure in the previous section. Fig. 6 shows the distribution of order parameter R s over different packings at P = 2 × 10 −2 . The evolution of the distribution with X S is completely different from that at P = 10 −3 . The distribution always has only a single peak at all X S , and the peak only continuously shifts with X S . In addition, Fig. 7 plots χ versus X S at the same pressure. Clearly, χ does not depend on the system size, which is a totally different observation in Fig. 3. These distinctive behaviors provide strong evidences for the lack of the first-order transition at higher pressure P = 2×10 −2 . Thus, there is R S X S =0.05 X S =0.08 X S =0.10 X S =0.12 Figure 6. Probability distributions of order parameter RS at P = 2 × 10 −2 . The distribution has only a single peak, and the peak continuously moves with XS. a critical point at a finite pressure P c between P = 10 −3 and P = 2 × 10 −2 . We repeated the above analysis at different pressures, plotted χ versus X S and observed the peak. At pressure P , we denote the peak height by χ peak (P ) and its position in X S -axis by X peak S (P ). In Fig. 8, χ peak (P ) versus P is plotted for various system sizes. χ peak (P ) depends on the system size at low pressure and not at high pressure. We determine that the first-order transition occurs at pressure P if χ peak (P ) exhibits the strong system size dependence even in our largest system size. Then, the first-order transition point is estimated as (P, X peak S (P )), which are shown as open symbols in Fig. 9. We also observe that χ peak (P ) increases with decreasing P in the high-pressure regime, which suggests that χ peak diverges at critical pressure P c . To precisely determine P c , we must simulate much larger systems, which is beyond the  scope of this work and can be addressed in the future. Here, we provide a rough estimate of P c as the highest pressure at which we observe the strong size effects in χ peak : P c ≈ 2 × 10 −3 , which is shown as a closed symbol in Fig. 9. In this phase diagram, the LS phase indicates the state with R S ≈ 1, where both components are jammed, and the L phase is the state with R S ≈ 0, where only the large component is jammed. At jamming point P → 0, the first-order transition is located at X S = 0.22, and the transition line continues for a finite pressure. This transition line is terminated at the critical point, which is estimated to be at X S = 0.195 and P = 2 × 10 −3 . line, and it terminates at the critical pressure. In this section, we explore the impact of this phase behavior on the mechanical properties of the system. We investigate bulk modulus B and shear modulus G of the system. Both moduli are calculated from the linear response formalism. Details of this calculation and formulation are provided in Appendix A. We present the results for N L = 64 in this section. Fig. 10 shows the values of the bulk and shear moduli averaged over different packing samples. Both moduli increase when X S increases. This result is consistent with the behaviors of R S sample in Fig. 5; the small particles participate in the rigid network and contribute to the rigidity of the system when X S increases. As in the case of R S sample , both moduli more dramatically increase when the pressure is decreased.
Similar to the analysis in the previous sections, we calculate the probability distributions of the bulk and shear moduli over different packings samples. Fig. 11 and Fig. 12 show the behaviors of the distributions of both moduli at high and low pressures. At P = 10 −3 , the distributions of both moduli are bimodal as in the case of R S . The peak at smaller moduli corresponds to the L phase, and the other corresponds to the LS phase. Thus, the elastic moduli also show discontinuous changes near the first-order transition in the jammed phase. Compare to the bulk moduli, the shear moduli have a broader distribution presumably because the fluctuation of the shear moduli becomes huge near jamming point P → 0; the bimodal distributions in shear moduli are smeared by the critical fluctuations from the jamming critical phenomena. At a higher pressure than the critical pressure, the bimodal distributions of both moduli are no longer observed. Therefore, the first-order phase transition is also characterized by the transition in mechanical properties. Finally, we discuss the pressure dependence of the moduli. In the monodisperse systems, many physical quantities are known to follow the critical power-law near the jamming transition. It has been established that the bulk modulus is independent of the pressure of the system, while shear modulus depends on the square roots of the pressure. However, it is not clear whether this jamming scaling is valid in the binary mixtures because of the discussed discontinuities of elastic moduli . In Fig. 13, we plot the average bulk modulus B sample and shear modulus G sample versus P with fixed X S . Clearly, there are two branches in the plot: one branch is associated with small X S , and the other is associated with large X S . At intermediate X S , we observe the jump from one branch to the other. In the lower branch, the systems are in the L phase, and the upper branch corresponds to the LS phase. In each branch, bulk moduli B are constant with pressure P , and shear moduli G depend on the square root of pressure P Hence, despite the discontinuities of the elastic moduli, the jamming scaling G ∝ p 1 2 holds in each jammed phase.

IV. CONCLUSION
In summary, we have numerically studied the structural and mechanical properties of the packings of binary mixture of particles with large size dispersity. We have presented strong evidences that the system exhibits the first-order phase transition between two distinct jammed phases: L phase, where only large particles are jammed, and LS phase, where both types of particles are jammed. This study was achieved by analyzing the statistics of the fraction R S of small particles that participate in the connecting network. The mean value of R S shows sharp increases at finite X S , and the increase becomes progres- sively steeper when the system size increases. The probability distribution of R S exhibits typical behaviors of the first-order transition, and the susceptibility associated to R S linearly increases with the system size. We have also shown that the elastic moduli can be used as the order parameter of this transition. The LS phase is more rigid than the L phase. The probability distribution of the bulk and shear moduli behave similarly to those of R S . Finally, we have shown that these moduli in each phase follow the critical power-law as in the monodisperse system. For nearly monodisperse systems, it has been established that the vibrational and transport properties exhibit the critical behavior near the jamming transition [3,4,6,7,12,[14][15][16]. As we established that the binary mixtures with large size dispersity exhibited another transition, it is interesting to study the vibrational properties of this system. We are now working in this direction.

ACKNOWLEDGMENTS
We thank K. Hukushima and H. Ikeda for useful discussions. We also thank B. P. Tighe for useful infor-mation on the result in two dimension This work was supported by JSPS KAKENHI Grants No. 18H05225, 19H01812, 19K14670, 20H01868, 20H00128. The computations were partially performed using the Research Center for Computational Science, Okazaki, Japan.

Appendix A: Formulation of elastic constants
In this appendix, we introduce the linear response formalism to calculate the elastic constants. We basically follow the formalism of Lemaitre [42]. U is a function of both particle positions r i and strain tensor η αβ . The stress tensor t αβ is defined as where D Dη αβ is the derivatives that impose the mechanical equilibrium f i = − ∂U ∂ ri = 0, and the second term vanishes because of this condition. The mechanical equilibrium leads to the following equation where ∂ 2 U ∂ rj ∂ ri = H ij is a 3 × 3 matrix with ij components of dynamical matrix H, and ∂ 2 U ∂η αβ ∂ ri = Ξ i,αβ is the nonaffine force field.
Elastic constants are defined as the second-order derivatives of the energy by the strain tensor η αβ and have following form This is decomposed into two terms as follows (A5) Unlike the stress tensor, the second term is non-zero under mechanical equilibrium and is called the non-affine correction to the elastic constants. The derivatives under constraint is replaced by the non-affine force field Ξ i,αβ , which leads to .
(A6) The second term is rewritten by diagonalizing dynam- Here, the nth eigen values and eigen vectors of H are represented as λ n and Ψ n , and the inner product of Ξ αβ and Ψ n is Ξ n αβ .