Cooling after shearing: three possible fates for dense granular materials

We perform discrete element simulations of freely cooling, dense granular materials, previously sheared at a constant rate. Particles are identical, frictional spheres interacting via linear springs and dashpots and the solid volume fraction is constant and equal to 60% during both shearing and cooling. We measure the average and the distributions of contacts per particle and the anisotropy of the contact network. We observe that the granular material, at the beginning of cooling, can be shear-jammed, fragile or unjammed. The initial state determines the subsequent evolution of the dense assembly into either an anisotropic solid, an isotropic or an anisotropic fluid, respectively. While anisotropic solids and isotropic fluids rapidly reach an apparent final steady configuration, the microstructure continues to evolve for anisotropic fluids. We explain this with the presence of vortices in the flow field that counteract the randomizing and structure-annihilating effect of collisions. We notice, in accordance with previous findings, that the initial fraction of mechanically stable particles permits to distinguish between shear-jammed, fragile or unjammed states and, therefore, determine beforehand the fate of the freely evolving granular materials. We also find that the fraction of mechanically stable particles is in a one-to-one relation with the average number of contacts per particle. The latter is, therefore, a variable that must be incorporated in continuum models of granular materials, even in the case of unjammed states, where it was widely accepted that the solid volume fraction was sufficient to describe the geometry of the system.

Grains in a silo are jammed until they start to flow like a gas in the proximity of the outlet [43]. Transitions from jamming to flow and from flow to jamming characterize the onset and the arrest of landslides and debris flows [13]. Also, in landslides or sediment transport in general, the solid particles typically flow over a substrate composed of the same flowing particles that can be mobilized, recently indicated to be in a fragile state [7]. Finally, it has been suggested that transitions from jammed/crystalline to fragile states can explain certain intermittent behaviour of intruders in granular materials [3,51].
Many works have been devoted to identify the most suitable state variables to characterize phase transitions in granular matter. A famous phase diagram for the transition between jammed and unjammed states [26], recently extended to include also shear-jammed and fragile states [9], makes use of the solid volume fraction and the shear stress. These variables have the advantage of being easily measured in experiments as well as in numerical simulations. On the other hand, their critical values at the transition cannot be derived from first principles, but need to be experimentally fit.
Indeed, rather than the solid volume fraction, the average number of contacts per particle, the coordination number, seems more fundamental, given its link to the number of degrees of freedom that need to be saturated to obtain a mechanically stable assembly of particles [2]. Unfortunately, even for the coordination number, its critical value at jamming is known a priori only for frictionless and infinitely frictional spheres: numerical simulations are needed for the intermediate cases [39]. It has also been suggested [2,9,22] that the fraction of non-rattlers grains, that is the number of particles with a minimum number of force-bearing contacts over the total number of particles, actually controls the transition between unjammed, fragile and shear-jammed states. It is intriguing that both the coordination number and the fraction of non-rattlers can be predicted, at least in the case of isotropic compression and shearing without rotation of frictionless particles, through a model that employs an history-dependent value of the critical volume fraction at jamming [22]. This seems to indicate that critical coordination number and critical solid volume fraction are, indeed, interchangeable concepts.
The role of the coordination number has been investigated also in a recent work [46], in which we have performed numerical simulations of unsteady, dense, homogeneous shearing at constant volume fraction. During the transient, the coordination number increased from zero to its value at steady state, and the granular material experienced the transition from an unjammed (fluid-like) to a shear-jammed (solid-like) behaviour at constant solid volume fraction. The phase transition was suggested to be characterized by a critical value of the coordination number, independent of the imposed volume fraction and the shear rate. The critical value of coordination number was found to be about the same that marks the development of rate-independent components of the stresses in steady, homogeneous shearing [12,41,46].
In our previous work [46], we distinguished the fluid-like from the solid-like behaviour based on the influence of the shear rate on the fluctuations of the coordination number during the transient. One question that we asked ourselves was whether the supposedly solid states retain their rigidity once we remove the applied shear and let the system relax.
Here, we perform discrete numerical simulations in which we, first, homogeneously shear an initially marginally jammed assembly of identical spheres. Then, we stop shearing at different times and let the system relax and dissipate its fluctuation kinetic energy (cooling). The value of the solid volume fraction is kept constant and equal to 0.6 during both shearing and cooling. The relaxation process of granular materials has been largely studied in the literature, and different volume fractions and initial conditions have been considered [27,30,32,36,37]. In this work, relaxation is allowed at different time steps of a time-evolving shearing simulation, so that the initial state of each cooling is characterized by the same volume fraction and shear rate, but different internal structure of the system. Depending on when we stop the shearing and start the cooling, we will show three different possible evolutions of the contact network and the velocity field, that we will associate to the dense granular material being, at the beginning of cooling, in a shear-jammed, fragile or unjammed state. We will also show that either the coordination number or the fraction of non-rattler grains are appropriate state variables to characterize the initial state and, therefore, determine the fate of dense granular materials.

Simulations
We use a freely available code, Mercury-DPM (www.mercu rydpm .org) [42,48], to perform discrete element simulations in a cubic box with 5000 identical spheres of mass density p and diameter d. The side of the box, L, is adjusted so that the solid volume fraction , i.e., the total volume of the particles inside the box divided by the volume of the box itself, has the desired value. The particles interact through linear springs and dashpots in both the normal and the tangential direction with respect to the line connecting the centres of the spheres in contact. The normal, e n , and tangential, e t , coefficients of collisional restitution are the negative ratios of post-to pre-collisional relative velocities of the interacting particles. Sliding occurs whenever the ratio of the normal to the tangential force at contact exceeds the friction coefficient, . The stiffness of the tangential spring, k t , is 2/7 the stiffness of the normal spring, k n . In all simulations, we employ e n = e t = 0.7 and = 0.3 . In the following, all variables are made dimensionless using the particle mass density, diameter and normal stiffness, so that velocities and times are given in units of k n d∕ p 1∕2 and p d∕k n 1∕2 , respectively. The simulations are performed using a dimensionless time step equal to 0.032, corresponding to 1/50th of the contact duration. However, data are stored at intervals of 1000 time steps, that is 32 units of dimensionless time, to reduce memory usage. We employ periodic boundary conditions and, after a preparation, perform a two-stage simulation procedure (shearing and cooling). In has been shown that a system size of 2000 particles (and here we employ 5000 particles) is sufficient to get size-independent measurements when periodic domains are employed [10,41]. We checked and confirmed the size-independence by performing the preparation and some of the subsequent simulations using both 2000 and 5000 particles. The simulations performed with 2000 particles (not shown here for brevity) are only characterized by a large scattering in the measured quantities.

Preparation
We prepare the granular system by following a widelyused three-step procedure [20,23], as illustrated in Fig. 1.
(P1) At first, we randomly place the spheres in the box at a solid volume fraction = 0.58 . In this stage, we use frictionless ( = 0 ) and elastic ( e n = 1 ) spheres, whose non dissipative nature guarantees homogeneous conditions. Then (P2), we set = 0.3 , e n = 0.7 and isotropically compress the box to increase the solid volume fraction to 0.6 and (P3) let the system to relax to a configuration where the residual kinetic energy of the particles is below a certain threshold. As explained in previous works [40,47], the rate at which the isotropic compression is performed determines the state of the assembly at the end of stage (P3).
Quick and slow compressions lead to jammed and unjammed states, respectively, despite the fact that the solid volume fraction is the same (0.6) and larger than the value at jamming (0.596 for friction 0.3, through interpolation of the measurements of Silbert [39] on granular packings). Unjammed states are characterized by zero coordination number, Z, and zero pressure, while jammed states by Z > 4 and nonzero pressure.

Shearing and cooling
After the preparation, the simulation procedure starts at t = 0 . In the first stage, we shear the system (S) by setting the upper and lower boundaries in motion with constant tangential velocities equal to V = 2.6 ⋅ 10 4 and −V , respectively, by employing Lees-Edwards [25] periodic boundary conditions along the direction perpendicular to the flow. Except for a short interval at the beginning of the stage (see later), the shearing is homogeneous and the dimensionless shear rate, ̇ , is equal to 3.16 ⋅ 10 −5 . We checked that decreasing the value of V (or that of ̇ ) by one order of magnitude, that is, equivalently, using stiffer particles, has no qualitative influence on the results. Starting from an unjammed state, at a constant volume fraction of 0.6, the granular material experiences, during shearing, a fluid-to-solid transition [4]. Starting from a marginally jammed state, instead, two phase transitions can be observed: a solid-to-fluid followed by a fluid-to-solid [47]. Here, we focus on the latter situation, because it provides a larger variety, in terms of properties of the contact network, of initial configurations for the last stage of our simulations.
Finally (C), at different times, t C , during the shearing, we stop the boundaries and let the system cool down, i.e., decrease its granular temperature-one-third of the mean square of the particle velocity fluctuations [14]-through inelastic collisions. Operatively, we take different configurations during the shearing as corresponding initial states for the cooling. We perform 34 cooling stages, whose characteristics are summarized in Table 1.
The sequence of stages employed in our numerical simulations, together with the associated imposed values of solid volume fraction and boundary velocity V, are depicted in Fig. 1.
We checked the homogeneity of the system during both shearing and cooling. The temporal evolution of the normalized profiles along the gradient direction, z, of the velocity in the flow direction, x, v x ∕V , and the local solid volume fraction, loc ∕0.6 , during shearing are plotted in Fig. 2a and b, respectively. The same profiles during the cooling simulations C34 and C14 are illustrated in Fig. 2c, d and e, f, respectively. Here and in what follows, the measurements are the results of coarse-graining [49].
During shearing (Fig. 2a), the instantaneous application of a relative velocity between the boundaries is not immediately transmitted within the system and the velocity profiles are initially nonlinear. After a short time interval ( t < 960 ), however, the shearing can be considered homogeneous [4]. Due to some difficulties of averaging in the proximity of the Lees-Edwards boundaries, measurements of volume fraction in regions of thickness 0.3L at the top and the bottom of the cubic box are not considered. In the rest of the domain ( 0.3 < z∕L < 0.7 ), the local solid volume fraction slightly varies along z, with deviations from its mean value 0.6 smaller than 1% (Fig. 2b).
During cooling, once the boundaries are arrested at t = t C , the motion of the particles suddenly stops in some cases (Fig. 2c), or slowly decelerates giving rise to a pronounced nonlinear distribution of the flow velocity in others (Fig. 2e). As explained in more details later, we use these behaviours to distinguish between initially jammed and unjammed states, respectively. In any case, the local volume fraction remains always rather homogeneous, with deviations smaller than 1% from 0.6 ( Fig. 2d, f).
We characterize the contact network through the fabric tensor, that gives a macroscopic measure of the contact orientation distribution. Sun and Sundaresan [41] defined the fabric tensor as a symmetric traceless tensor: where I is the identity matrix, N tc is the total number of force-carrying contacts and ( ) is the unit contact normal vector of the th contact. Here, we prefer to use the definition of the fabric tensor, , suggested by Weinhart et al. [49] and valid in the case of mono-dispersed particles, which has a nonzero trace equal to the product of the solid volume fraction and the coordination number. In our configuration, is simply related to through During shearing, as already noticed [41], the only non-negligible off-diagonal component of the fabric tensor is F xz , where, as mentioned, x and z are the flow and gradient directions, respectively (y is the vorticity direction). This component represents, therefore, a good measure of the degree of anisotropy of the contact network [15,20].

Results
We report the temporal evolution of the coordination number and the shear component of the fabric tensor during shearing (stage S) in Fig. 3 (black crosses). As previously stated, we take t = 0 at the beginning of shearing, when the system is isotropic ( F xz = 0 ) and jammed (Z slightly less than 5, but well above 4.5, i.e., roughly the minimum value at jamming obtained from interpolating the data of Silbert [39]). As soon as the shearing is imposed, the coordination number decreases and the fabric anisotropy increases. Around t = 2000 , there is a sharp decrease in the coordination number, which reaches a minimum of 0.5. This first (solidto-fluid) transition is accompanied by a sharp decrease in the fabric anisotropy. Then, both Z and F xz increase over time and, eventually, when the coordination number is large enough, the granular material manifests a rate-independent behaviour as shown by Vescovi et al. [46], a signature of a second (fluid-to-solid) transition. More details about the above mentioned transient during shearing and a proposal for a mathematical continuum model to reproduce it are reported elsewhere [47].
As mentioned, here we formally interrupt the shearing at different times t C during the transient and let the system cool down (stage C). The states at the beginning of cooling are marked as symbols of large size in Fig. 3, whereas small symbols refer to the subsequent evolution during cooling. in the x-z plane, respectively, for the three different behaviours (rows). Columns from left to right show the evolution of the configurations from the beginning of cooling, t = t C , to the end of simulation, t = t F (after 16,000 units of time, i.e. t F = t C + 16, 000 ). We identify three different behaviours of the system during cooling. For the sake of clarity, for each kind of behaviour, we plot only a couple of representative cooling simulations in Fig. 3 and only one in Figs. 4 and 5. The same considerations apply nonetheless to all the simulations reported in Table 1. At the largest values of the initial   coordination number (red circles in Fig. 3a), irrespective of the anisotropy of the fabric (Fig. 3b), the contact network quickly reaches a steady state in which both the coordination number, always larger than 4 (signature of solid configuration), and the fabric anisotropy are constant and slightly less than the values at the beginning of cooling (Fig. 3). This behaviour is associated with initially Shear-Jammed (SJ) states in which a strong contact network that spans the entire domain persists during cooling (Fig. 4). The presence of this contact network causes the almost immediate arrest of all the particles once the shearing is interrupted (Fig. 5). For a narrow range of initial coordination numbers (blue triangles in Fig. 3a), no matter what the anisotropy is (Fig. 3b), the granular material experiences a catastrophic collapse of the entire contact network and rapidly reaches a seemingly steady, isotropic state with almost zero coordination number (signature of fluid configuration), the homogeneous cooling state of granular gases [14]. At first, the contact network spans the entire domain (Fig. 4), as for initially Shear-Jammed states, and is strong enough to almost immediately arrest the motion of the particles in the x-direction after shearing is removed (Fig. 5). The contact network is, however, Fragile (F) and therefore unstable to perturbations, such as those associated with the particle velocity fluctuations that still persist during cooling. After the collapse, most of the contacts are wiped out (Fig. 4).
In some cases, the state at the beginning of cooling is Marginally Shear-Jammed (MSJ), a critical condition between Shear-Jammed and Fragile (purple diamonds in Fig. 3). The contact network first relax to an apparently solid configuration, that is, however, unstable. Indeed, after a certain amount of time, the network suddenly collapses as for initially Fragile states (Fig. 3).
At the lowest values of the initial coordination number (green squares in Fig. 3a), once again with no apparent dependence on F xz (Fig. 3b), the system quickly reaches and UnJammed states (simulations C34, C31 and C14, respectively). The plots are shaded according to kinetic energy, from low (light) to high (dark) non negligible, but less than unity, values of the coordination number (if the coordination number at the beginning of cooling, Z C , was not already in that range) and small, but nonzero, values of the fabric anisotropy. Then, both Z and F xz slowly decays with time, without apparently reaching a steady state (Fig. 3). We were not able to observe a steady configuration of the contact network even when we increased of an order of magnitude the total simulation time (not shown here for brevity). We associate the above described behaviour with initially UnJammed (UJ) states in which the contact network is much weaker than in Shear-Jammed or Fragile states (Fig. 4). Indeed, the particle velocity persists even after the removal of shearing (Figs. 2e and 5). Due to the inelastic, randomizing collisions, the contact network progressively disappears (Fig. 4), accompanied by the development of long-lasting vortices (Fig. 5) (sometimes, present since the beginning of cooling), as previously observed in freely cooling granular gases at low and high volume fraction [27,30] and inside shear bands [44]. The presence of these meso-structures slows down the energy dissipation and the obliteration of the contact network, although a final isotropic state at vanishingly small coordination number will likely be asymptotically reached. Table 1 summarizes the values of the coordination number at the beginning ( t = t C ) and at the end ( t = t F ) of the cooling stage, Z C and Z F , respectively, for the entire set of simulations. In the next section we measure the distributions of contacts during cooling and propose a quantitative criterion to classify the different initial states. representative cases of initially Shear-Jammed (C34), Fragile (C31), Marginally Shear-Jammed (C33) and UnJammed (C14) states. Although the quantitative characterization of the contacts, as discussed in what follows, certainly depends on the properties of the particles (friction and restitution coefficients) and their size distribution, we do expect that the qualitative trends are universal. The distribution of contacts is measured through the contact fraction:

Distribution of contacts during cooling
where N = 5000 is the total number of particles and N c is the number of particles having c contacts, for c = 0 and 1. Also shown is the case c < 2 , for which the contact fraction f c coincides with f R , the fraction of rattlers, i.e., particles having less than 2 contacts which makes them mechanically unstable [1]: with The fraction of rattlers is trivially related to the fraction of non-rattlers, which, as already mentioned, was previously suggested to control the transition between unjammed, fragile and shearjammed states in 2D configurations [2,9]. We point out that the choice of the maximum number of contacts in the rattler definition is somehow arbitrary for particles with finite friction. Other possibilities, such as including also particles with 2 and 3 contacts among rattlers, would not qualitatively alter the following considerations.
It seems reasonable that systems with a larger fraction of non-rattlers, or, equivalently, a smaller fraction of rattlers, behave more like a solid than like a fluid and are more stable against perturbations. Figure 6a indicates that, for initially Shear-Jammed states, the initial fraction of rattlers is lower than roughly 5% (the initial fraction of nonrattlers is higher than roughly 95%), and the fractions of contacts remain nearly constant during cooling. Almost all rattlers are floaters (particles with zero contacts), whereas the contribution of particles with one contact is negligible.
Instead, for initially Fragile states (Fig. 6b), the starting fraction of rattlers is about 5%, but almost immediately jumps to 100%, mainly consisting of floaters, in accordance with the annihilation of the contact network shown in Fig. 4. There is a certain amount of particles having one contact which, however, slowly decreases over time and reaches a seemingly steady value of 0.2%.
For initially Marginally Shear-Jammed states (Fig. 6c), the starting fraction of rattlers is about 5% (the initial fraction of non-rattlers is about 95%) and slightly increases during cooling until it reaches a value of about 10%. Then, there is a sudden jump to 100% at t − t C ≈ 2500 , like for initially Fragile states. Before and after the transition, the qualitative distributions of the fractions of contacts resemble those in initially Shear-Jammed (Fig. 6a) and Fragile (Fig. 6b) states, respectively. The fact that f R immediately before the jump in the number of floaters in Fig. 6c is higher that the initial fraction of rattlers in initially fragile states (Fig. 6b) suggests that there is not a unique critical value of f R , or, equivalently, f NR , at which there is a solid-to-fluid transition. We conjecture that the stability of the contact network is also influenced by the amount of fluctuation kinetic energy in the system (a measure of it is the granular temperature). More energetic systems need larger fractions of non-rattlers to remain stable (or, equivalently, remain stable at smaller fractions of rattlers). Given that the granular temperature decreases over time during cooling, this might explain why granular materials destabilize at larger values of f R in the case of initially Marginally Shear-Jammed systems.
For initially UnJammed states (Fig. 6d), the fraction of rattlers is higher than 75% and slowly increases during cooling up to 90%. Unlike the case of Fig. 6b, there is a significant amount of particles with one contact besides floaters. Also, it is clear that the system is still evolving and has not yet reached a steady configuration of contacts.
The values of fraction of non-rattlers at the beginning and at the end of cooling, f C NR and f F NR , for the entire set of simulations, are reported in Table 1. As in previous works [2,9], we found a one-to-one relation between the coordination number and the fraction of non-rattlers during both shearing and cooling (Fig. 7). This is in accordance with the suggestion that a common hidden state variable determines both Z and f NR [22] and is also somehow reassuring: the fraction of non-rattlers has a clear physical link with the stability of granular materials in a solid-like state, but the coordination number is a variable that is probably easier to incorporate in continuum models that include the microstructure evolution of the contact network [41]. In Fig. 7, we also show that the initially fragile states cover a narrow range of values of both Z and f NR . However, as mentioned, the granular temperature of the system at the beginning of cooling is likely to have an influence on that range.

Conclusions
We have performed discrete element simulations of identical, soft, frictional spheres in which an isotropically compressed cubic box with a fixed number of particles was first sheared at a constant rate. Then, the shear was suddenly removed and the granular material progressively dissipated its kinetic energy. Despite the fact that the solid volume fraction was kept constant and equal to 0.6, we have observed that the granular material, during cooling, can evolve into either (i) an anisotropic solid, (ii) an isotropic or (iii) an anisotropic fluid, depending on the characteristics of the contact network once shearing is interrupted.
Anisotropic solids, characterized by coordination numbers higher than four and nonzero shear components of the fabric tensor, derive from initially Shear-Jammed (SJ) states, in which the fraction of mechanically stable particles exceeds 95%. The chains of particles in contact span the entire domain and induce the rapid arrest of the flow after the removal of shearing. Their persistence explains the solidlike behaviour of the assembly.
Isotropic fluids, in which the coordination number and the anisotropy of the fabric tensor are practically zero, originate from initially Fragile (F) states, in which the fraction of mechanically stable particles is in the range 92-94%. Although the contact network spans the entire domain and is able to rapidly arrest the particle flow during cooling, the residual kinetic energy associated with the particle velocity fluctuations is seemingly sufficient to destabilize the network and cause its catastrophic collapse. The subsequent, freely evolving state strongly resembles the homogeneous cooling state of granular gases, although, to our knowledge, the latter have never been observed at such large volume fractions.
Anisotropic fluids, in which the coordination number is less than four and the anisotropy of the fabric is nonzero, follow initially UnJammed (UJ) states, in which the fraction of mechanically stable particles is less than 90%. In this case, the contact network cannot rapidly arrest the flow during cooling. We have observed that vortex instabilities develop and are able to slow down the energy dissipation and the destruction of the contact chains by random, inelastic collisions. As a consequence, we have not been able to reach a final, steady configuration for these anisotropic fluids.
We have determined a one-to-one relation between the coordination number and the fraction of mechanically stable particles. Hence, the thresholds for the classification of the initial states, which the fate of granular materials during cooling depends on, can be rephrased in terms of average number of contacts per particle, a more commonly adopted variable in continuum models based on the microstructure evolution of the fabric.
The present work reinforces the idea that the solid volume fraction is not enough to determine the state of granular materials and the microstructure of the contact network must be included in mathematical models, at least when dealing with unsteady problems. Although this might have been obvious to some in situations where the particle stresses are rate-independent [11,19,24,41], it has been recently shown that it is crucial also near the fluid-solid transition, where the particle stresses have also a rate-dependent component [22,46]. Here, we have shown that even when the granular material behaviour is far from solid-like, the volume fraction is not sufficient to characterize the geometry of the assembly, at least in dense situations.
It remains to investigate, in more detail, (i) the influence of the granular temperature on the range of existence of initially fragile states and, more generally, on the different granular dynamics during cooling; and (ii) the roles played by the solid volume fraction and the particle properties (coefficients of sliding friction and collisional restitution).
In any case, reproducing the results of the present simulations will be a severe test for future, comprehensive models of granular materials. provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.