A computational study of transient shear banding in soft jammed solids

TL;DR

This study uses 3D numerical simulations of soft spheres to investigate transient shear banding during the yielding of jammed soft solids, revealing robust features unaffected by damping mechanisms or boundary conditions.

Contribution

The paper introduces a detailed computational analysis of transient shear banding in jammed soft solids, highlighting the robustness of stress-overshoot phenomena across different simulation setups.

Findings

Stress-overshoot and transient banding are robust features in overdamped systems.

Boundary conditions and damping mechanisms influence load curves and velocity profiles.

Transient shear banding occurs consistently regardless of simulation parameters.

Abstract

We have designed 3D numerical simulations of a soft spheres model, with size polidispersity and in athermal conditions, to study the transient shear banding that occurs during yielding of jammed soft solids. We analyze the effects of different types of drag coefficients used in the simulations and compare the results obtained using Lees-Edwards periodic boundary conditions with the case in which the same model solid is confined between two walls. The specific damping mechanism and the different boundary conditions indeed modify the load curves and the velocity profiles in the transient regime. Nevertheless, we find that the presence of a stress-overshoot and of a related transient banding phenomenon for large enough samples are a robust feature for overdamped systems, where their presence do not depend on the specific drag used and on the different boundary conditions.

Equations

$$G^{\prime }(\omega ,{\gamma }_0)=\frac{\omega }{{\gamma }_0\pi }{\int }_{t_0}^{t_0+2\pi /\omega }{\sigma }_{xy}(t)sin(\omega t)dt,$$

$$G^{\prime \prime }(\omega ,{\gamma }_0)=\frac{\omega }{{\gamma }_0\pi }{\int }_{t_0}^{t_0+2\pi /\omega }{\sigma }_{xy}(t)cos(\omega t)dt$$

$$m\frac{d^2{\vec{r}}_i}{d{t}^2}=-{\zeta }_{DPD}\sum _{j(\ne i)}\omega (r_{ij})({\hat{r}}_{ij}\cdot {\vec{v}}_{ij}){\hat{r}}_{ij}-{▽}_{{\vec{r}}_i}U$$

$$m\frac{d^2{\vec{r}}_i}{d{t}^2}=-{\zeta }_{SD}\left(\frac{d{\vec{r}}_i}{dt}-\dot{\gamma }{y}_i{\vec{e}}_x\right)-{▽}_{{\vec{r}}_i}U$$

Full text

I Introduction

Solids whose microscopic constituents are densely packed into an amorphous assembly and basically insensitive to thermal fluctuations form an integral part of our everyday life, with examples ranging from pharmaceutical, cosmetic and food products to wet cement. Hence understanding the flow properties of such materials is of technological importance Bibette et al. (2003); Van Hecke (2009); Bonn et al. (2017). Under shear deformations, these materials flow only when the applied shear stress is above a threshold value, the yield stress. As the imposed strain increases at a fixed shear rates, they show an initial linear increase in the stress, often followed by a stress overshoot, beyond which the system yields before reaching a steady flow state Varnik et al. (2004); Fielding (2014); Divoux et al. (2016); Shrivastav et al. (2016). The steady state shear stress ($\sigma$) as a function of the applied shear rate ($\dot{\gamma}$) provides a constitutive behavior (flow curve) that in most cases is well described by the Herschel-Bulkley (HB) curve $\sigma=\sigma_{Y}+\kappa\dot{\gamma}^{n}$, where the yield stress $\sigma_{Y}$, the effective viscosity $\kappa$ and the exponent $n$ are in principle material specific Hébraud and Lequeux (1998); Nicolas et al. (2014). While the constitutive behavior under flow is interesting in itself and is the subject of many investigations Goyon et al. (2008); Lemaître and Caroli (2009); Tighe et al. (2010); Martens et al. (2012); Seth et al. (2011); Vasisht et al. (2018a); Nicolas et al. (2018), here we are interested instead in the transient behavior of the material before it reaches the steady-state flow, where it maintains features of both its solid and fluid response. Before reaching the steady flowing state, jammed soft solids often display flow instabilities which manifest in terms of strong spatial inhomogeneities in the flow profile, often called shear localization or shear banding, even when the material is homogeneously driven.

Shear banding in complex fluids like polymer solutions and wormlike micelles has been extensively studied and, in many cases, its origin is understood in terms of flow alignment of the microstructure or flow induced crystallization, through theoretical approaches that couple the flow fields with the microstructure Dhont (1999); Olmsted (2008); Dhont and Briels (2008); Lerouge and Berret (2009); Bonn et al. (2017). In jammed and soft glassy materials, these phenomena are the subject of intense investigations Bécu et al. (2006); Fielding et al. (2009); Divoux et al. (2010); Besseling et al. (2010); Mansard et al. (2011); Martens et al. (2012); Fielding (2014); Irani et al. (2014); Divoux et al. (2016); Bonn et al. (2017); Gross and Varnik (2018) and their possible association to an underlying non-equilibrium phase transition is increasingly debated Moorcroft et al. (2011); Adams et al. (2011); Wisitsorasak and Wolynes (2017); Parisi et al. (2017); Ozawa et al. (2018); Popović et al. (2018). However a good understanding of their microscopic origin is fundamentally lacking and particularly elusive when the shear banding is observed as a transient feature that does not persist in steady-state Fielding (2014); Divoux et al. (2010, 2016); Shrivastav et al. (2016).

Microscopic computer simulations can be effective in bringing new insight, but most of the existing studies have been performed in 2D and in the quasi-static limit (i.e. zero shear rate) while experiments are always performed at finite shear rates and nearly always in 3D. Together with meso-scale simulations based on elasto-plastic models Martens et al. (2012); Nicolas et al. (2018), most of existing microscopic simulation studies have focused on the part of the phenomenon that can be rationalized in terms of the emergence of local plastic rearrangements in an elastic background Lemaître and Caroli (2009). That is, they do not capture the coupling of the microscopic dynamics to the imposed deformation rate. Recent works, instead, have pointed precisely to the crucial role of the rate dependent dynamics, calling for dedicated numerical investigations Divoux et al. (2010, 2011); Fielding (2014).

To be able to address such questions, we have devised a numerical study using 3D microscopic simulations and finite deformation rates. In this paper we analyze the outcomes for the shear start-up response under varying shearing protocols and boundary conditions, with the goal to extract the robust features that do not depend on the specifics of the numerical simulations. We study a jammed suspension (at volume fraction $\phi\simeq 0.7$) of polydisperse spheres with soft repulsive interactions, i.e., a deeply jammed soft amorphous solid. We find that the shear start-up is characterized by a stress-overshoot and formation of a transient band, which eventually disappears upon reaching the steady state. We extensively test whether such findings depend on different implementations of the viscous drag due to the background solvent or on different types of boundary conditions. Our study indicates that the stress-overshoot and the shear banding during the shear start-up are robust features emerging from the microscopic physics of our model material, they do not qualitatively change with the different shearing protocols used here and are mainly controlled by the sample age that we can change through the preparation protocol.

The paper is organized as follows. In section II we discuss the model system and the sample preparation protocol as well as the rheological and structural characterization of the initial sample configurations. In section III we discuss the shearing protocols followed in our work for performing finite shear rate simulations. Sections IV and V contain the results comparing the transient rheological properties obtained from different protocols. At the end, in section VI we present conclusions and discussion.

II Model and sample preparation

The 3D model for soft amorphous material we investigate is composed of a non-Brownian suspension of soft repulsive particles interacting via a truncated and shifted Lennard-Jones potential Weeks et al. (1971), given by $U(r)=4\epsilon\left[(a_{ij}/r_{ij})^{12}-(a_{ij}/r_{ij})^{6}\right]+\epsilon$, for $r_{ij}\leq 2^{1/6}a_{ij}$, else $U(r_{ij})=0$. Here $\epsilon$ is the unit energy in the simulations, $r_{ij}$ being the center to center distance between the particle $i$ and $j$ and $a_{ij}=0.5(a_{i}+a_{j})$, with $a_{i}$ and $a_{j}$ being the diameter of particles $i$ and $j$ respectively. The diameters of the particles are drawn from a Gaussian distribution with variance of 10%, whose mean is used as unit length $a$. All the simulations are performed at a volume fraction $\phi\approx 70\%$, consisting of $10^{5}$ (97556) particles (unless otherwise mentioned). Albeit simple, this type of model has been successfully used for numerical simulations of soft solids in various contexts, and proven to capture the fundamental physics of their behavior under deformation Varnik et al. (2003); Irani et al. (2014); Shrivastav et al. (2016); Vasisht et al. (2018a).

We prepare the numerical samples to be sheared using the following procedure. An initial FCC crystal at the chosen volume fraction of $0.70$ (with $lx=ly=lz=42.1798a$) is melted at $T=5.0\epsilon/k_{B}$ and equilibrated at the same temperature using NVT Molecular Dynamics (MD) for $\simeq 5\cdot 10^{4}$ MD steps, with a MD timestep of $\Delta t=0.001\tau_{0}$ (where $\tau_{0}=\sqrt{ma^{2}/\epsilon}$ is the unit of time and $m$ is the unit mass). After making sure that there is no signature of crystallinity in the equilibrated sample by measuring the bond orientational order parameter $Q_{6}$ Steinhardt et al. (1983), the melt is subjected to a systematic temperature quench. From the initial temperature of $T=5.0\epsilon/k_{B}$ we decrease the system temperature by $\Delta T$, after which we let the system relax at the temperature $T-\Delta T$ for $10^{4}$ MD steps. We repeat this procedure several times, until $T=0.001\epsilon/k_{B}$ is reached. By changing the $\Delta T$ we control the cooling rate $\Gamma$. The samples discussed here have been prepared for $\Gamma$ corresponding to $5\cdot 10^{-2},5\cdot 10^{-3},5\cdot 10^{-4}\epsilon/(k_{B}\tau_{0})$. After the system reaches $T=0.001\epsilon/k_{B}$, we perform an energy minimization using the conjugate gradient (CG) method to take the system to the zero temperature limit. For the fastest quench rate case, we minimise the energy directly from the equilibrium liquid state at $T=5.0\epsilon/k_{B}$ using CG. Following the above described protocol we prepare five independent samples at each cooling rate (except for the lowest one, where we have three samples). The initial configurations so obtained are subjected to shear deformation at finite shear rate. Before studying the shear start-up behaviour we carefully analyze various properties of the initial sample configurations.

In Fig. 1(a) we show the pressure measured as a function of the temperature as the sample is quenched at a chosen cooling rate $\Gamma$. The pressure dependence on the temperature does not change with $\Gamma$ until $T=3.0\epsilon/k_{B}$, which in our model corresponds to the onset of slow glassy dynamics. Below such onset temperature, the pressure achieved through our cooling protocol decreases with decreasing $\Gamma$. In Fig. 1(b), we plot the energy of the closest local minimum accessible to the system, obtained through CG, as a function of $T$. One can clearly see that, upon decreasing $\Gamma$, we access deeper local minima (or inherent structures) of the total potential energy Sastry et al. (1998); Stillinger (1995); Mosayebi et al. (2014); Ashwin et al. (2004) below the onset temperature. We next characterize the mechanical and structural properties of these inherent structure configurations obtained from difference $\Gamma$.

II.1 Linear viscoelastic response of the initial configurations

In order to characterize the mechanical properties of the initial configurations, we compute the complex modulus by performing small amplitude oscillatory shear. By applying a shear strain $\gamma(t)=\gamma_{0}sin(\omega t)$, for a strain amplitude of $\gamma_{0}=1\%$, we monitor the stress response for varying frequencies $\omega$ Larson (1999). We monitor the energy and pressure evolution during the oscillatory shear cycles and once the system reaches a saturation in these quantities as a function of the number of cycles, we extract the linear visco-elastic moduli (respectively storage and loss) using

[TABLE]

In Fig. 2 we show the $G^{\prime}(\omega,\gamma_{0}=0.01)$ and $G^{\prime\prime}(\omega,\gamma_{0}=0.01)$ as a function of the applied frequency $\omega$ for the initial configuration prepared at a cooling rate $\Gamma=5\cdot 10^{-4}\epsilon/k_{B}\tau_{0}$. The data show that the sample is solid, with $G^{\prime}>>G"$ at low frequencies and the zero frequency storage modulus can be extracted by extrapolating the data to $\omega=0$. We compute the zero frequency storage modulus for samples prepared with different $\Gamma$, as shown in the inset of Fig. 2. The arrow in the inset of Fig. 2 corresponds to the fastest cooling rate, where a liquid configuration is directly subjected to CG minimization. With decreasing the cooling rate $\Gamma$, the configurations in deeper local minima of the potential energy landscape, as identified in Fig. 1 (b), clearly correspond to solids with higher mechanical strength.

II.2 Structural analysis

In order to gain insight into the microscopic structural underpinnings, we analyze the structure of the samples by constructing a 3D radical Voronoi tesselation with the Voro++ open source software library Rycroft (2009) and compute the Voronoi polyhedra. This method is well suited for our samples made of polydisperse spheres. For samples prepared at different cooling rate $\Gamma$, we obtain the statistics of various Voronoi polyhedra (defined by number of faces, edges and vertices). In Fig. 3, we show the fraction of different polyhedra (which constitute more than $5\%$ of each sample) as a function of the cooling rate. The data show that the fraction of Voronoi dodecahedra is clearly affected by the cooling rate $\Gamma$, whereas the contribution of the other polyhedra does not change much. A Voronoi dodecahedron corresponds to a local packing which is an icosahedron, where neighbouring particles form an icosahedron around a reference particle. Relatively higher percentages of local icosahedral packing are expected and widely observed in slowly quenched supercooled liquids Steinhardt et al. (1981), as well as in glasses, with similar spherically symmetric interactions Mosayebi et al. (2012); Royall and Williams (2015); Ronceray and Harrowell (2015).

Having characterized the rheological and structural properties of the initial configurations prepared at the various cooling rates, we now analyze the rheological response of a sample obtained with the slowest cooling rate, by subjecting it to a continuous shear deformation at a finite rate.

III Shearing protocols

The samples characterized as described above are subjected to a shear deformation of strain amplitude $\gamma(t)$ with constant imposed rate $\dot{\gamma}$, either using Lees-Edwards boundary conditions (LEBC) or using a wall-based protocol (WB). Note that in all simulations discussed below [x,y,z] refers respectively to flow, gradient and vorticity directions. In the WB protocol, we confine the sample in $y$ direction by freezing two layers of particles at the simulation box boundary. $L_{y}$ is the width of the box in the $y$ direction for the confined samples. All shear rates are measured in units of $\tau_{0}^{-1}$.

For LEBC we solve the equations of motion in two different ways. In the first way, which in the following we call LEBC1, we use the following dissipative particles dynamics (DPD) equations of motion Nicolas et al. (2016):

[TABLE]

where $m$ is the mass of the particle and the first term in the right hand side (RHS) is the damping force which depends on the damping coefficient $\zeta_{DPD}$. The relative velocity $\vec{v}_{ij}=\vec{v}_{j}-\vec{v}_{i}$ is computed over a cut-off distance $r_{ij}\leq 2.5a_{ij}$, with the weight factor $\omega(r_{ij})=1$. These choices for the cut-off distance and the weight factor are consistent with other studies in the literature for similar systems Xu et al. (2005); Maloney and Robbins (2008); Puosi et al. (2015); Nicolas et al. (2016). In the viscous damping we have only considered the radial contribution of the relative velocities, since the particles are point-like and the main sources of change in the velocity are the inter-particle forces which are purely radial Puosi et al. (2014). Nevertheless, the tangential contribution to the damping forces can also be included as $((\hat{r}_{ij}\cdot\vec{v}_{ij})\hat{r}_{ij})+(\vec{v}_{ij}-(\hat{r}_{ij}\cdot\vec{v}_{ij})\hat{r}_{ij})$ Baumgarten and Tighe (2017). In most of this work we have used the formulation with only the radial contribution, but in one of the following sections we briefly discuss the effect of including the transverse term, which has also been investigated recently in Irani et al. (2019). The second term in the RHS is the force due to the interactions between particles.

In the second way, which we indicate as LEBC2, we solve equations of motion where the solvent drag is Stokes-like and only depends on the absolute particle velocity Ikeda et al. (2013); Colombo and Del Gado (2014); Hoh and Zia (2015); Su et al. (2017), as in a free draining approximation, given by

[TABLE]

where again $m$ is the particle mass and the first term in the RHS is the damping force which depends on the damping coefficient $\zeta_{SD}$. Hence the drag force is proportional to the difference between the particle’s velocity $d\vec{r}_{i}/dt$ and an affine background velocity, dictated by imposed shear rate, given by $\dot{\gamma}y_{i}\vec{e}_{x}$. The second term in the RHS is the force due to the interactions between particles. While one concern with this second approach is that it does not strictly conserve momentum (i.e., it is not Galileian-invariant) Maloney and Robbins (2008); Soddemann et al. (2003), its use can still be justified in systems where most of the stress induced through the imposed deformation is due to inter-particle interactions and the contribution of the solvent is a minor correction, since this type of drag term allows for simpler and faster simulations, with no need to adjust additional parameters such as the cut-off or the form of the weight factor in the DPD one. We will show in the following that this free draining approximation works well in the case of deeply jammed systems as those of interest here.

For the wall based shear deformation tests, which we call WB in the following, we follow the procedure of Varnik and co-workers in Varnik et al. (2004). Two walls confine the samples along the direction $\hat{y}$, at a relative distance $L_{y}$: one wall moves at a velocity $\vec{v}=v_{x}^{wall}\hat{x}=\dot{\gamma}L_{y}\hat{x}$), while the other is kept fixed. The particles that form the wall are completely frozen during the evolution of the system, but the interactions between the wall and the sample particles are the same as in the sample. For the WB shear deformation we use only the DPD approach as in eqn.(3). In fact, if the drag coefficient is directly proportional to the particle velocity as in eqn.(4), in the WB simulations the particles near the moving wall feel more drag than the ones away from it. As a consequence, for large samples the time needed for particles near to the non-moving wall to sense the deformation is quite a long (i.e., quite larger than the sound speed in the system). This behavior could be mistaken for a shear banding but it is only a numerical artifact due to the use of eqn.(4) in this geometry and disappears using instead the DPD equation, eqn.(3).

Finally, both eqs. 3 and 4 include the particle inertia since this allows us to use the same efficient and precise algorithms devised for MD Frenkel and Smit (2001). Nevertheless we can effectively study the dynamics in the overdamping limit regime by suitably adjusting the damping coefficients $\zeta_{SD}$ and $\zeta_{DPD}$. As a measure of the extent of damping, we can define an inertial quality factor $Q=\tau_{\mathrm{d}amp}/\tau_{\mathrm{v}ib}$, that measures the ratio between the time scale over which the inertial motion is damped $\tau_{\mathrm{d}amp}=m/\zeta$, due to the solvent drag, and the time scale $\tau_{vib}=\sqrt{ma^{2}/\epsilon}$ over which a particle of mass $m$ and diameter $a$ is accelerated by the force $\epsilon/a$. For a fully over-damped system $Q\rightarrow 0$. It has been shown that for the athermal conditions considered here, and for similar type of interactions, $Q\approx 1$ well approximates the overdamped regime Nicolas et al. (2016); Vasisht et al. (2018b). Hence here we focus on $Q\approx 1$ and show also that the results we obtained do not vary significantly by decreasing $Q$ or varying it around $1$. As stated above, the shear rates in the following are always expressed in units $\tau_{0}^{-1}\equiv\tau_{\mathrm{vib}}^{-1}$.

All simulations have been performed using LAMMPS molecular dynamics package Plimpton (1995), with modifications to the source code to incorporate continuous polydisperse interactions, since we use a distribution of particle sizes and the interaction shape depends on the particle diameters. In the Fig. 4, we show a rendering of samples prepared at $\Gamma=5\cdot 10^{-4}\epsilon/(k_{B}\tau_{0})$ and sheared at $\dot{\gamma}=10^{-4}\tau_{0}^{-1}$ (the snapshots are both taken at $\gamma=0.12$, during the stress decay) using LEBC and wall based protocols. We use the color scheme to show the velocities of the particles which vary between 0 $a/\tau_{0}$ (white) and 0.005 $a/\tau_{0}$ (black).

In the following section we present a comparative study of the rheological response starting from the load curves, using the two different boundary conditions and the two different equations of motion introduced above.

IV Load curves and flow profiles: comparison of the different shearing protocols.

We begin with the rheological data obtained for a sample prepared at cooling rate $\Gamma=5\cdot 10^{-4}\epsilon/(k_{B}\tau_{0})$, i.e., a well annealed sample, and for a sample obtained by directly minimizing the energy of the high temperature fluid with CG, i.e., a poorly annealed sample. Both samples are sheared at $\dot{\gamma}=10^{-4}\tau_{0}^{-1}$ and $\dot{\gamma}=10^{-2}\tau_{0}^{-1}$, using the different protocols discussed above. All data in this section corresponds to the overdamping limit, with the inertial quality factor $Q=1$. Since all the simulations are performed in athermal conditions, the shear stress is computed from the virial stress tensor $\sigma_{\alpha\beta}=\frac{1}{V}\sum_{i}\sum_{j>i}r_{ij}^{\alpha}f_{ij}^{\beta}$, where V($=lx*ly*lz$) is the volume of the system, $r_{ij}$ represents the distance between particle $i$ and $j$, $f_{ij}$ is the force on the particle $i$ due to particle $j$ and ${\alpha,\beta}\in{X,Y,Z}$. We indicate the shear component of the stress $\sigma_{XY}$ with $\sigma$ and the virial pressure is obtained as $P=\frac{1}{3}\left(\sigma_{XX}+\sigma_{YY}+\sigma_{ZZ}\right)$. The first normal stress difference is computed as $\sigma_{11}=\sigma_{XX}-\sigma_{YY}$.

In Fig. 5 we show load curves, i.e., the shear stress $\sigma$ vs. applied strain $\gamma$, for the three different protocols (LEBC1, LEBC2 and WB) and for the two different samples (the well annealed one and the poorly annealed one mentioned above), sheared at $\dot{\gamma}=10^{-4}\tau_{0}^{-1}$, a relatively low shear rate, and $\dot{\gamma}=10^{-2}\tau_{0}^{-1}$, a higher shear rate. At low shear rate, all the protocols show comparable load curves. For the well annealed sample, independent of the protocol used, we find that the initial linear regime is followed by a stress overshoot before reaching a steady state value (see Fig. 5 (a)). The poorly annealed sample do not show any overshoot (see Fig. 5). At relatively high shear rates (see Fig. 5 (b) and (d)), load curves consistently show the dependence on preparation protocol w.r.t the stress overshoot, but we also note that the dependence on shear protocol become more prominent. We also note that, approaching the steady state, the stress fluctuations are larger in the case of LEBC2 compared to LEBC1 and WB. When comparing LEBC1 and WB at high rate, we note that the difference in the boundary conditions does not seem to affect the value of the stress overshoot nor of the steady-state value, which instead change with changing the damping mechanism (LEBC1 vs LEBC2).

Now we focus our discussion on the low shear rate case ($\dot{\gamma}=10^{-4}\tau_{0}^{-1}$) and on the well annealed samples, to analyze the corresponding velocity profiles. Figs. 6, 7 and 8 show the velocity profiles respectively obtained with the LEBC1, LEBC2 and WB protocols from the same initial configurations. For different regimes in the load curves (Figs. 6 (a), 7(a) and 8(a)), we show the associated velocity profiles in Figs. 6, 7 and 8 (b) (i-ix). The velocity profiles are computed over a strain window of $\Delta\gamma=0.02$. In spite of some differences in the details of the velocity profiles, we can recognize the following features for all the shearing protocols considered. 1) In the linear response regime of the load curve (i-ii), the velocity profile show a homogeneous flow with local shear rates (obtained as the slope of the velocity profile) similar to that of applied shear rate. 2) In the vicinity of the stress overshoot (iii) a flow instability develops, as indicated by the deviation from the linear velocity profile. Following the stress overshoot, the stress decays back, with $d\sigma/d\gamma<0$ and in this region we find a back-flow in the system as the velocity profile show a negative slope Moorcroft and Fielding (2013); Fielding (2014). We note that the width of the shear band seems to be set by the amount of back flow in the system (iv). 3) Upon further shearing, the width of the band grows and eventually the whole system flow homogeneously. The shear component of the stress has a much weaker (although non-negligible) dependence on the strain and the pressure evolution with $\gamma$ clearly shows that saturation of the stress tensor only happens once the profile returns to be homogeneous (Figs. 6, 7 and 8 (c)) Olmsted (2008); Divoux et al. (2010).

When analyzing the local packing through the Voronoi tessellation, we find that the fraction of icosahedrally packed particles has a similar trend. Fig. 6 (d) shows, for the LEBC1 protocol, how the fraction of icosahedrally packed particles evolves with the increasing strain and the time required for it to saturate is consistent with the vanishing of the transient shear banding. The connections between the evolution of the local packing and that of the shear inhomogeneities, as well as their dependence on the shear rate have been thoroughly investigated in Vasisht et al. (2017). We also find that the evolution of the system during the shear banding is accompanied by a positive first normal stress difference ($\sigma_{XX}-\sigma_{YY}$), indicating the presence of dilation, which, although quite noisy, also seems to grow as the bands develop and saturate as the flow becomes homogeneous (see Fig. 6 (e)). When comparing further the different protocols across Figs. 6 - 8, we note that in WB the shear bands always nucleate near the walls and that the time taken to reach an homogeneous flow state is slightly longer than in LEBC1 and LEBC2.

We conclude that the features identified above, the stress overshoot, the band formation and disappearing, its correlation with the evolution of the normal stresses and with the changes in the local packing, are robust across the different protocols utilized. Hence they must not be the result of numerical artifacts and are instead inherent of well aged jammed solids sheared at sufficiently low rate.

In the next section we discuss the possible dependence of our findings on the value of the drag coefficient chosen.

V Dependence on dissipation co-efficient

The contribution of inertial terms to the rheology of jammed suspensions has been extensively explored in steady-state in Salerno et al. (2012); Nicolas et al. (2016); Vasisht et al. (2018b). Here we address instead how the stress overshoot and the shear banding during the transient preceding the steady state is affected by the dissipation co-efficient. To do this, we consider the LEBC1 protocol and vary the value of the coefficient $\zeta_{DPD}$ between $2$ ($m\tau_{0}^{-1}$) and $0.01$ ($m\tau_{0}^{-1}$). This means that the inertial quality factor $Q$ defined above varies from 0.5 to 100, i.e. from the case where the damping time is much shorter than the characteristic time of the inertial motion (i.e., the overdamped limit) to a case where it is  100 times longer (underdamped case). In Nicolas et al. (2016); Vasisht et al. (2018b) it has been shown that a strongly underdamped dynamics may qualitatively change the flow properties of the material, with the possibility to develop even non-monotonic flow curves. This is not the case here, where the flow curves remain monotonic, suggesting that for the whole range of $Q$ considered here the high density and the jammed conditions control the material flow, rather than the particle inertia. In Fig. 9 (a) we show the load curves of the same well annealed sample studied with LEBC1, but now with different values of the inertial quality factor Q. The applied shear rate corresponds to $10^{-4}\tau_{0}^{-1}$. The data show that the initial linear response remains unaffected and the stress overshoot seems to slightly decrease with the increasing Q. The decay after the overshoot becomes increasingly steeper with increasing $Q$.

If we express the applied shear rate in terms of Weissenberg number (${\mathrm{W}_{i}}=\zeta_{DPD}a^{2}/\epsilon\dot{\gamma}$), the higher the quality factor, the lower is $W_{i}$. As a consequence, one could think of the increasing quality factor $Q$ as a way to reach effectively lower shear rates. One might expect a trivial scaling of the load curves if the load curve is presented in terms of $\sigma$ vs. $\gamma/\dot{\gamma}$, with $\dot{\gamma}$ in terms of Weissenberg number. But since the initial linear response regime is unaffected by the extent of damping, this scaling would not work. The effect of the inertial contribution can be more intricate, with inertia playing a role in increasing the kinetic temperature in the athermal system and hence leading to a softening of the system during shearing Nicolas et al. (2016); Vasisht et al. (2018b).

In Figs. 9 (b)-(e), we show the velocity profiles related to different $Q$ and computed at a strain $\gamma\approx 0.2$ by averaging over a strain window $\Delta\gamma=0.02$. In all the cases, as the stress decays from the overshoot, we observe the formation of shear bands.

Finally, in the data of Fig. 9 (a) we also show the results obtained in the protocol LEBC1 when including the transverse contribution to the DPD drag, as discussed in section III, to demonstrate that such modification changes slightly the value of the overshoot and the shape of the decay towards the steady state, but it does not modify the general phenomenology of formation of shear band observed (see also the profile shown in Fig. 9 (c)).

With the idea to explore further how the specific form of drag used may change the results obtained, we note that in LEBC1 another important variable corresponds to the pair-wise dissipation cutoff used in the DPD drag. For the results shown so far, we have chosen such cut-off to be $2.5\mathrm{a_{ij}}$, as done in Nicolas et al. (2016). In Fig. 10 (a), we plot together the load curves obtained for LEBC1 with the DPD cut-off $2.5a_{ij}$ and $1.5a_{ij}$, along with the results for the LEBC2 protocol. We observe that the linear response regime is unaffected and so is the values of the stress overshoot. The stress decay from the the overshoot shows instead a dependence on the cut-off chosen. We observe that the results obtained with the LEBC2 protocol approach those for LEBC1 if the DPD cut-off distance goes down towards a particle diameter, consistent with the fact that LEBC2 corresponds to the free draining approximation. The related velocity profiles are shown in Figs. 10 (b)-(d) (these are computed at $\gamma=0.2$ and averaged over $\Delta\gamma=0.02$), indicating that the decay of the overshoot is always associated to the formation of transient shear bands for both values of the DPD cut-off.

Finally, we have also analyzed how the difference or similarities just described depend on the shear rate. We summarize the outcome of this study in Fig. 11, which shows, for the well annealed sample, the difference between the stress overshoot and the steady state value of the shear stress as a function of the shear rate, for different protocols as well as for two different damping coefficients for LEBC1. The data indicate that changing $Q$ in the range of values explored here does not significantly affect the results at sufficiently low rates, whereas one should expect to see significant differences upon increasing the shear rate. Note that we have expressed the shear rate in terms of $\tau_{0}^{-1}$, instead of using a viscous time scale, in order to compare different protocols where viscous terms are handled in different ways. Different protocols show that the occurrence of the stress overshoot (which is essential for the formation of shear bands) is robust and its dependence on the steady-state stress across different protocols is similar in the whole range of shear rates (the lines through the data in the figure are power law fits). Hence the stress overshoot and the banding observed at low rates are robust features of the systems we are studying, they are genuinely representative of the emerging physics and not the result of artifacts in the numerical simulations.

VI Dependence on the sample age

At the end this study, we want to emphasize that the tendency to have flow inhomogeneity upon yielding is determined mainly by the age of the samples (which determines the stress overshoot and eventual decay to steady state), due to the different degree of frozen-in stresses. In the context of the study reported here, we have analyzed the dependence of the stress overshoot on the shear rate and of the sample age. Fig. 12 shows the difference between the stress overshoot and the steady state value of the shear stress as a function of the shear rate, for different cooling rates $\Gamma$, corresponding to different sample ages. The results clearly show how changing the sample age can qualitatively change the overshoot and its rate dependence. We have investigated how the age of the samples also determine the persistence of the flow inhomogeneities, not only their presence, in the companion letter Vasisht et al. (2017). We refer the reader to that paper for further insights.

VII Conclusions

We have devised a 3D numerical study of a jammed suspension of soft spheres, polydisperse in size, under shear. In particular here we have explored different choices for imposing the shear deformation and boundary conditions. We have compared the use of Lees-Edwards boundary conditions with simulations where the samples are confined within walls. We have also compared the use of a DPD drag term to the free draining approximation. Finally, we have compared simulations with different degrees of inertia, quantified through the inertia quality factor $Q$. The comparison has been done in terms of load curves and velocity profiles during the transient that leads to the steady-state flow. In all cases, we find that at low rates the shear stress develops an overshoot followed by a relatively long decay (not necessarily gradual) towards the steady state value. Such phenomenon is associated to a transient banding with a part of the material that is basically stuck, and the rest flowing. On the basis of the extended comparison carried on here, we propose that these features (the stress overshoot and the transient shear banding) are the genuine results of the emerging response of the material upon yielding and not the consequences of numerical artifacts or unphysical choices in the simulations parameters. We have been able to elucidate the dependence of the stress overshoot on the different choices for the simulation parameters and to clarify that what controls the response reported here is ultimately the age of the samples, which determines the spatial distribution of frozen-in stresses in the initial solid, as it has to be expected for jammed out-of-equilibrium materials.

Acknowledgements. The authors thank Georgetown University and the Swiss National Science Foundation (Grant No. PP00P2 150738) for support.