On the stress-overshoot in cluster crystals under shear

Using non-equilibrium molecular dynamics simulations we study the yielding behavior of a model cluster crystal formed by ultrasoft particles under shear. We investigate the evolution of stress as a function of strain for different shear rates, $\dot{\gamma}$, and temperatures. The stress-strain relation displays a pronounced maximum at the yielding point; the height of this maximum, $\sigma_{\rm p}$, increases via a power law with increasing shear rate and tends to saturate to a finite value in the limit shear rate goes to zero (at least within the considered temperature range). Interestingly, this behavior can be captured by the Herschel-Bulkley type model which, for a given temperature, allows us to predicts a static yield stress $\sigma^{0}_{\rm p}$ (in the shear rate tending to zero limit), a characteristic timescale $\tau_{c}$, and the exponent $\alpha$ of the power-law decay of the $\sigma_{\rm p}$ at high shear rates. Furthermore, for different temperatures, the $\sigma_{\rm p}$ can be scaled as functions of $\dot{\gamma}$ onto a single master curve when scaled by corresponding $\tau_{\rm c}$ and ${\sigma}_{\rm p}^{0}$. Moreover, for a given shear rate, $\sigma_{\rm p}$ displays a logarithmic dependence on temperature. Again, the $\sigma_{\rm p}$-$T$ curves for different shear rates can be scaled on a single logarithmic master curve when scaled by corresponding fitting parameters.

Using non-equilibrium molecular dynamics simulations we study the yielding behavior of a model cluster crystal formed by ultrasoft particles under shear. We investigate the evolution of stress as a function of strain for different shear rates,γ, and temperatures. The stress-strain relation displays a pronounced maximum at the yielding point; the height of this maximum, σp, increases via a power law with increasing shear range and tends to saturate to a finite value in the limit shear rate goes to zero (at least within the considered temperature range). Interestingly, this behavior can be captured by the Herschel-Bulkley type model which, for a given temperature, allows us to predicts a static yield stress σ 0 p (in the shear rate tending to zero limit), a characteristic timescale τc, and the exponent α of the power-law decay of the σp at high shear rates. Furthermore, for different temperatures, the σp can be scaled as functions ofγ onto a single master curve when scaled by corresponding τc and σ 0 p . Moreover, for a given shear rate, σp displays a logarithmic dependence on temperature. Again, the σp-T curves for different shear rates can be scaled on a single logarithmic master curve when scaled by corresponding fitting parameters.

I. INTRODUCTION
Cluster crystals [1] represent an archetypical ordered phase encountered in a particular class of soft matter systems, so-called ultrasoft particles, whose interaction attains at vanishing interparticle distance a finite energy penalty: such cluster phases can be found in systems if the Fourier transform of the potentials shows negative components [2]. At low and intermediate densities, a disordered phase of clusters of overlapping particles is found, characterized by a rather polydisperse distribution in the cluster occupancy. If at a fixed temperature the density is increased this disordered phase transforms via a first order phase transition into a cluster BCC phase and then -upon further increasing the density -into a cluster FCC phase. In these ordered particle configurations the lattice sites of the BCC and of the FCC lattice are occupied by clusters of overlapping particles which are now rather monodisperse in their cluster occupancy [1,[3][4][5][6][7]. Even though the fact that mutually repelling particles are able to form stable clusters is at least at first sight counter-intuitive, the physics behind this intriguing phase is meanwhile well understood in terms of classical density functional theory [8,9].
Later contributions in literature have been dedicated to the equilibrium dynamics of cluster crystals [7,10]. Again, these systems show an intriguingly new behaviour: the longtime dynamics is diffusive which is the result of hopping particles from one cluster to a neighboring cluster, the distribution of the jump length is found to be exponential at short distances. At large distances, this distribution follows a power-law decay, i.e., a behavior reminiscent of Lévy flights [11]. As expected, the diffusion coefficient exhibits the Arrhenius behavior con-firming that the hoping of particles is an activated process. This particular type of diffusion of particles makes cluster crystals an appropriate model system to study the mechanical response of defect-rich crystals [12], notably the effect of the defect-dynamics on the yielding of crystalline solids, an area which is less explored in the recently developed theories to understand the deformation of crystalline solids [13,14].
Only recently, out-of-equilibrium investigations have been dedicated to cluster crystals [10,15,16] . In a recent computer simulation study on the steady-state rheology of cluster crystals an intriguing shear-induced selfassembly, was observed, namely a shear-induced string formation is observed at very high shear rates. When the shear rate is further decreased string-like structures appear and at very low shear rates a shear shear-banding regime occurs.
This contribution is dedicated to the yielding behavior of cluster crystals, a feature which, so far, has not been investigated at all. In particular, we investigate the effect of shear rate and temperature on the mechanical failure of these materials. In this work, we study with the help of extensive computer simulations the yielding behavior of the cluster crystals under steady shear conditions by examining the response of the stress as a function of strain. Similar as most of the preceding investigations on cluster forming systems, our investigations are based on the socalled generalized exponential model interaction with index n (GEM-n), a potential which combines the bounded nature of the interaction (required to guarantee clustering) with a functionally simple mathematical form that is easily amenable to computer simulations. We will show that the stress exhibits a maximum after an initial linear increase; the height of this maximum depends on the shear rate and temperature.
The paper is organized as follows: in Section II, we introduce the model and details of the simulations and the related protocols. Results are presented and discussed in Section III, while the final Section contains conclud-ing and summarizing remarks and an outlook to future, related investigations.

II. SIMULATION DETAILS
In our cluster crystal system particles interact via the so-called generalized exponential (GEM-n) potential [1], assuming in this contribution n = 4; this interaction is defined via where d and set the length-and energy-scales of the model, respectively; note that we use -in contrast to the usual notation in previous, related contributionsthe symbol d for the range of the interparticle interaction, since the conventional symbol σ is reserved in this manuscript to denote the stress.
In the simulations the GEM-4 potential is truncated and shifted to zero at a distance r c = 2.2d. The units of temperature (T ), density (ρ), and time (t), are given by k B T / , ρd 3 , and t 0 = d m/ , respectively. Here, m is the mass of particles and k B is Boltzmann's constant. In our simulations we set the values of , d, m, and k B equal to unity.
All investigations have been carried out with the LAMMPS simulation package [17]: we have performed molecular dynamics (MD) simulations in an NVT ensemble of N = 3328 particles (with periodic boundary conditions), integrating Newton's equations-of-motions with the Verlet velocity integration scheme [18]; a time increment of ∆t = 0.005t 0 has been used. The temperature was kept constant with the help of a DPD thermostat which preserves hydrodynamics [19]. For each state point 50 independent simulation runs have been performed; observables were then obtained by averaging the related data over these runs.
The system has been studied at a fixed density, namely ρ = 6.5, and at four different temperatures, i.e. T = 0.8, 0.75, 0.7, 0.65, and 0.6. At these state points the system is in a stable FCC cluster phase, where each site of the FCC lattice is occupied by a cluster of overlapping particles (see phase diagram shown in [1]). According to the results available in literature [1,5,11] the average number of particles N c pertaining to a cluster ranges for the considered state points around the value N c 13.
The initial configurations for our simulations are ideal FCC cluster crystals where each lattice site is occupied by N c = 13 completely overlapping particles and assuming a lattice distance that is compatible with the given value of ρ. Starting from this configuration the system is equilibrated independently at each of the above specified temperature values over 10 6 MD steps. A snapshot of the resulting equilibrium configuration (at ρ = 6.5 and T = 0.6) is shown in panel (a) of Fig. 1.
We shear the bulk cluster crystal by imposing a planar Couette flow in the (x, z)-plane and use Lees-Edwards boundary conditions [21]. To this end a constant shear rate is applied along the x-direction; thus, in this setup the z-and y-directions are the gradient-and vorticitydirections, respectively. The shear rate,γ, considered in our investigations the simulations ranges from,γ = 10 −5 toγ = 10 −1 (given in units of t −1 0 ). A schematic illustration of the setup of our shearing experiment is presented in panel (b) of Fig. 1: here the front view of the snapshot (i.e., its projection onto the (x, z)-plane) of the equilibrated cluster crystal (as displayed in panel (a) of Fig. 1) is shown; the velocity gradient is indicated by the blue arrows.
We note that in contrast to Ref. [15], we consider here rather low shear rates for the following reasons: (i) the DPD thermostat (as implemented in the LAMMPs package) becomes in our investigations unstable at shear rates larger thanγ = 10 −1 ; (ii) the present work focuses on an understanding of the yielding behavior of cluster crystals in the limit of vanishing shear rates tends to zero; still it should be mentioned that, on the other hand, also very low shear rates (i.e.,γ 10 −5 ) are difficult to access, as related simulations become computationally expensive.

III. RESULTS
In an effort to understand the yielding of the cluster crystal under shear, we record the evolution of the stress, σ xz (t), during the shearing experiment as a function of strain,γt, assuming different values for the shear rateγ.
From the quantities available during the simulation run the stress σ xz (t) is calculated via the Irving-Kirkwood expression [22]: Here, V represents the total volume of the system, m is the mass of the particles, v i,x (t) and v i,z (t) represent the x-and z-components of the velocity of particle i, r ij,x (t) is the x-component of the displacement vector between particles i and j, and F ij,z (t) denotes the zcomponent of the force between particles i and j. The angular brackets in relation (2) represent an averaging over the aforementioned 50 independent runs. In the following subsections we discuss the impact of the shear rate and of the temperature on the stress.
A. Impact of the shear rate on the stress In panels (a) and (b) of Fig. 2 we display the timeevolution of the stress, σ xz (t) , as a function of the strainγt for five different shear rates, namelyγ = 10 −1 , 10 −2 , 10 −3 , 10 −4 , and 10 −5 at temperatures T = 0.6 (panel (a)) and T = 0.8 (panel (b)); similar graphs also exist for the other temperatures investigated, but are not on display here. For all values ofγ we observe that the stress first increases with the shear rate and then reaches a pronounced maximum. At this peak in the stress-strain response curve the cluster crystal yields and the stress decays beyond this maximum. We observe that this decay becomes sharper (eventually even abrupt) as the shear rate is decreased. Beyond this decay essentially two different archetypical scenarios of the stress-strain curve can be identified: (i) for moderate shear rates (i.e., forγ 10 −1 ) the curve levels off without any significant characteristic features; (ii) for smaller shear rates, however, the stress-shear curves show characteristic secondary peaks, which become more and more pronounced as the shear rates decreases. These additional peaks can be attributed to partial release events of the stress via local particle arrangements, a feature which becomes more pronounced as the shear rate is lowered. When after such a local rearrangement event the strain increases, the cluster crystal yields again and another secondary peaks appears in the stress-strain curve which has a lower height than the preceding (and the primary) peak.
The data presented in panels (a) and (b) of Fig. 2 also provide evidence that at a fixed temperature, the height of the above mentioned peak in the stress, denoted henceforward as σ p , decreases as the shear rate is lowered. This also indicates that the resistance applied by the cluster crystal to the shear forces reduces as the shear rate is lowered. This observation is visualized in panel (c) of Fig. 2, where we have plotted σ p as a function of the shear rateγ for all the temperatures investigated. In an effort to quantify this effect and to provide thus more insight into the observed phenomena we have fitted σ p (γ) via the following functional form, known in literature as the Herschel-Bulkley type expression [23,24], which was originally developed for the steady state, but which we extend now to the out-of-equilibrium case: The maximum in the stress-strain curve σ p is also referred to in literature as the "static yield stress" [24]. In case of amorphous solids, it has be observed in computer simulations that the stress overshoot shows a logarithmic dependence on the shear rate [25,26]. By fitting the simulation data to the above expression one can extract a value for σ 0 p , i.e., the value of the static yield stress asγ tends to zero; we conclude that σ 0 p will thus correspond to the yield stress which is defined as the minimum stress required to initiate plastic deformation in the cluster crystal [25]. We note that alternatively, this quantity could also be calculated from constant stress simulations [25], an issue which we postpone to future investigations. The temperature-dependent parameter A(T ) in the above equation is related to a characteristic timescale, τ c (T ), which we will introduce below in Eq. (4). Finally, α is the exponent that characterizes the power-law decay of the stress overshoot as a function oḟ γ.
The solid lines in panel (c) of Fig. 2 show the simulation data, along with the curves which emerge by fitting these results by the Herschel-Bulkley type expression defined in Eq. (3). The emerging values of σ 0 p (T ) and A(T ) are listed in Table I. The exponent α is found to be essentially independent of the temperature in the considered range of temperatures, namely α 0.43. It is worth mentioning that a related exponent α occurs also in the Herschel-Bulkley model for the flow curve (i.e., the variation of the steady-state stress as a function of shear rate). We leave investigations related to this issue to later contributions.   an essentially linear decay with increasing temperature, a feature that can probably be attributed to the fact that for a cluster crystal at equilibrium the number of particles that hop from one cluster to a neighbouring one decreases with decreasing temperature. Further, the different curves σ p = σ p (γ), as shown in panel (c) of Fig. 2 can be mapped onto one single master curve by scaling the shear rate by a temperaturedependent timescale, τ c (T ), the latter one being defined as τ c = A/σ 0 p 1/α ; the values of this timescale for different temperatures are accumulated in the Table I. Starting from Equ. (3) this master curve has thus the form [27]: which the different σ p (γ)-curves can be mapped. Such a scaling behaviour provides evidence of a universal scenario of yielding in cluster crystals which -at least in the considered range of temperatures -is independent of the temperature.

B. The effect of temperature
It can be seen from panel (c) of Fig. 2 that the yielding of cluster crystals depends strongly on the temperature. In the following we analyze this dependence by investigating the behavior of σ p as a function of temperature for different shear rates. Data shown in panel (a) of Fig. 4 provides evidence that σ p decreases logarithmically as the temperature increases for all the shear rates considered. Thus the data can be fitted via where B (γ) and C (γ) are suitably chosen parameters. The dashed lines in panel (a) of Fig. 4 represent the fitting function defined in the above equation. The related coefficients, B (γ) and C (γ), shown in Fig. 4(b), are themselves functions of the shear rate: (i) B (γ) decreases at high shear rates withγ, it saturates at low shear rates; (ii) theγ-dependence of C (γ) does not show any significant features. Finally we point out that the quantity [B(γ) − σ p (T )]/D(γ) = ln T can be mapped onto a single master curve, shown in panel (c) of Fig. 4 as a black dashed line. It should be noted that such a logarithmic decay of σ p as a function of temperature is also observed in thermal glasses [25,26,28]. This is the more astonishing as such systems are metastable, thus the age of the sample becomes an important feature. Cluster crystals, on the other hand, are equilibrium systems, and aging in this case is not a relevant phenomenon.

IV. SUMMARY AND OUTLOOK
For this contribution we have made an extensive simulation-based study to understand the impact of the shear rate and of the temperature on the yielding of cluster crystals which is exposed to shear forces. In particular we investigated the behavior of the stress overshoot for different shear rates and temperatures by analysing the stress-strain curve. We have found that for a given temperature the height of the stress overshoot, σ p , increases in a power-law manner with increasing shear rate in a Herschel-Bulkley type manner; the related exponent is essentially independent of temperature -at least for the range of the temperatures considered in this work. At low shear rates, σ 0 p tend to saturate to a finite value, σ 0 p . However, we note here that a systematic finite-size analysis is needed to ensure the non-vanishing value of σ 0 p .
Further, the behavior of the stress overshoot as a function of shear rate is found to have a universal nature, thus different curves of σ p (γ) can be mapped onto a single master curve. The static yield stress is found to decay in a linear manner with temperature. At a given shear rate, the height of the stress overshoot exhibits a logarithmic dependence on the temperature. This behavior is also found to be universal, thus different σ p (T ) curves can be mapped onto a single master curve. This logarithmic behavior is a characteristic feature of defect-rich systems such as thermal glasses [25,26] and is motivated by the Ree-Eyring viscosity theory [23,29].
V. ACKNOWLEDGMENTS