Role of thermal expansion heterogeneity in the cryogenic rejuvenation of metallic glasses

Cryogenic rejuvenation in metallic glasses reported in Ketov et al 's experiment (Nature(2015)524,200) has attracted much attention, both in experiments and numerical studies. The atomic mechanism of rejuvenation has been conjectured to be related to the heterogeneity of the glassy state, but the quantitative evidence is still elusive. Here we use molecular dynamics simulations of a model metallic glass to investigate the heterogeneity in the local thermal expansion. We then combine the resulting spatial distribution of thermal expansion with a continuum mechanics calculation to infer the internal stresses caused by a thermal cycle. Comparing the internal stress with the local yield stress, we prove that the heterogeneity in thermomechanical response has the potential to trigger local shear transformations, and therefore to induce rejuvenation during a cryogenic thermal cycling.


I. INTRODUCTION
Glasses are solid materials, disordered at small scale, but homogeneous and isotropic at large scales. It has been understood for some time, however, that their elastic properties are not uniform if they are computed on mesoscale. This heterogeneous elasticity has been characterized by a number of simulations of model glasses, more or less realistic, and also by a variety of experimental tools [1][2][3][4]. The heterogeneity in elastic properties plays a key role for understanding vibrational spectra, thermal transport, but also far from equilibrium behavior such as plastic deformation. Macroscopic, plastic deformation proceeds via microscopic rearrangements, the "shear transformation", that take place preferentially on "weak spots" with typically low elastic properties. A large plastic deformation, if it can be performed homogeneously, tends to "rejuvenate" the system, by bringing it to states of higher internal energy.
Recently, it was discovered that an alternative to this mechanical rejuvenation could be obtained, in some metallic glasses (MGs), by thermal cycling towards low temperatures [5]. This is quite surprising, as from a purely thermal point of view one would expect that bringing the system to a lower temperature does not populate higher energy states, and should only result in slow aging. One possible interpretation of the effect, suggested by Hufnagel [6], is that the rejuvenation is due to the creation of internal stresses during the temperature cycling. The origin of such stresses could be twofold: local heterogeneities in the temperature field due to the very fast heating process which were investigated by some of us in a recent molecular dynamics study [7]. Alternatively, an heterogeneous thermal expansion coefficient, similar to the heterogeneity in elastic constants, would also create internal stresses even if the temperature field is homogeneous. If the corresponding mechanical stresses are large enough, they would generate plastic activity internal to the material, resulting in rejuvenation. Molecular dynamics simulations, performed on relatively small samples and using very fast cycling rates combined with artificial thermalisation methods, do not allow us to precisely quantify the origin of the internal stresses and rejuvenation induced by the cycling process.
In this work we take an alternative route that bridges the gap between molecular dynamics and experiments, by using the molecular dynamics simulations to obtain quasi equilibrium properties and assessing the importance of the induced mechanical stresses by using a con-tinuum description. Our first step is to quantify the heterogeneity in thermal expansion.
We therefore present the first quantitative study of the local thermal expansion in a model metallic glass, putting the thermoelastic properties of glasses on the same footing as their elastic properties. We use the resulting spatial distribution of the thermal expansion as an input in a continuum model that allows us to estimate the magnitude of the stresses generated by differential thermal expansion. These stresses are finally compared to the local yield stresses that must be overcome in order to trigger plastic activity. Our main conclusion is that, while the differential thermal expansion is indeed capable of generating plastic activity and thermal rejuvenation upon cryogenic cycling, the effect is unlikely to be observable at the scale of molecular dynamics simulations.

II. THE LOCAL THERMAL EXPANSION COEFFICIENT
The thermal expansion coefficient quantifies the volume strain caused by a change in temperature, and can be written for a macroscopic system as where v is the volume strain induced by the temperature change ∆T . In metallic glasses, its order of magnitude is around 10 −5 K −1 [8].
In order to define a local thermal expansion coefficient, we use the coarse graining approach originally proposed by Goldhirsch and Goldenberg [9,10], that allows one to define continuum fields in the sense of hydrodynamics or elasticity, starting from atomic positions.
The same approach was used previously to define local elastic constants, with results consistent with those obtained with different methods. The approach starts with the definition of a meso scale displacement field u(r, t) [10] from the atomic positions: where u i (t) is the displacement of atom i at time t, starting from a reference position, and φ(x) is the coarse graining function, in this paper we choose a Gaussian φ(x) with σ the coarse graining scale and v ef f = (4πσ 2 ) 3/2 .
The local strain is calculated, under the assumption of small deformations, as ij (r) = . We obtain the local thermal expansion coefficient (LTEC) from the local volumetric strain as where α(r) is the LTEC at position r, and v (r) is local volumetric strain caused by the temperature change ∆T .
We apply this procedure to a model metallic glass Cu 50 Zr 50 described by an embeddedatom method(EAM) potential [11]. The molecular dynamics(MD) simulations were conducted by an open source classical MD software: LAMMPS [12]. The simulation samples contained 8000 atoms, and metallic glass samples were obtained by quenching from 2000 K to 1 K with two different quenching rates, 10 9 K/s and 10 13 K/s, respectively. We collected configurations at 1 K , and then reheated the samples to 100 K Configurations at 1 K and 100 K were used to calculate the local thermal expansion coefficients from the coarsegrained volumetric strain. During all the process, the system was controlled in pressure and temperature (NPT ensemble) using a Nose-Hoover thermostat [13] and a Parrinello-Rahman Barostat [14], which was always maintained at 0 bar. As shown in Fig.1(a),(b), the probabil- ity distribution function of the LTEC is typically Gaussian, with a standard deviation that increases as the coarse graining size decreases, which is qualitatively consistent with similar calculations of the local elastic moduli [2,9]. As the coarse grain size decreases, the thermal response, similar to the mechanical response, becomes heterogeneous [9]. Interestingly, Fig.1(c) shows that the macroscopic thermal expansion coefficient is only weakly dependent on thermal history, but its probability distribution function is much more sensitive to the quenching rate. High quenching rates result in higher energy states, with more heterogeneity in thermal properties.
In the following we will be particularly interested in discussing the possibility that the stress field induced by the heterogeneous thermal expansion induces shear transformations (ST), i.e. localized irreversible plastic events that are recognized as the elementary building blocks of plasticity in amorphous materials [15][16][17]. Previous work on similar systems indicates that the shear transformations involve a few tens of atoms [18,19], hence we will concentrate on results obtained for a coarse graining size σ = 5Å, with a coarse graining volume that contains around 30 atoms. Results obtained for different coarse graining sizes, and the influence of the value of σ, are also investigated in the discussion section.

III. INTERNAL STRESS BY THERMAL EXPANSION HETEROGENEITY
When the metallic glasses are heated, due to the heterogeneity of LTEC, internal stresses are generated by the mismatch in the local thermal deformation. Extracting these thermal stresses directly from molecular dynamics simulations is not easy, in particular as the system can also evolve during thermal cycling due to local, thermally activated processes that have no relation to the differential thermal expansion. In order to estimate the role of the latter, we will adopt an alternative approach based on the assumption that the material remains a linear elastic one and can be described, at the coarse grained level, using continuum theory, as previously demonstrated in experiments [3] and simulations [9]. In the following we will use the simplifying assumption that the only heterogeneity is the one of the thermal expansion, while elastic properties (bulk and shear modulus) are uniform. While this is clearly a simplification, we could not detect any correlation between heterogeneity in local elastic constants and in thermal expansion. We also note that, on the scale under consideration, the variance of the local elastic moduli is small compared to the mean (see Fig.A1 in Appendix), so that we expect this assumption to provide at least a correct order of magnitude of the local stresses.
With these assumptions, the effective external stress generated by the local thermal expansion is K∆T α(r)δ ij where K is the bulk modulus, ∆T is the temperature change and α(r) is the LTEC at location r. δ ij the Kronecker symbol. At mechanical equilibrium, the elastic stress would balance the thermal stress σ ij as using small deformation conditions and linear elasticity, the equations for the displacement field are where u i is the i component of the displacement field at location r. Since the thermal expansion external force is a conservative vector field, we search the solution for u i in the can be simplified as Integrating in a finite system with periodic boundary conditions, equation 6 can be written whereᾱ ≡ α(r)dV V . Equation 7 can be solved numerically using as an input the LTEC obtained from the MD simulations, and the resulting elastic stress can be calculated from the potential field Φ(r). We characterize the resulting stress tensor using the von Mises stress [20] where σ i are the eigenvalues of the tensor. As shown in Fig. 2, the probability distribution function p(σ th ) shows a long tail distribution which we have fitted by the expression p L (x) = (γ/η)(1 + x/η) −γ−1 (Lomax distribution [21]). We do not have any particular justification for this fit, which is only used for convenience in the following.

IV. PROBABILITY OF TRIGGERING AN ATOMIC REARRANGEMENT
If the internal stress is larger than the local yield stress, we expect that thermal cycling will trigger plastic activity. In order to estimate the corresponding density of events, we first need to estimate the distribution of the local yield stresses in the material. We employ GPa;B=108.92,107.25 GPa for 10 9 K/s and 10 13 K/s, respectively which is also shown in Appendix and ∆T = 100 K the method described by Patinet and coworkers [22] sometimes described as "frozen matrix" method [2,22,23]. In a low temperature system, a sphere of finite radius r c is deformed in different directions while imposing a purely affine deformation of the surrounding system, until a local plastic instability is observed, in the form of a jump in the stress-strain curve.
The von Mises stress just before the first jump in the stress strain curve is used to define a local yield stress. The coarse graining function implicit in this approach is a discontinuous one, which separates the affinely deforming medium from the rearranging sphere. In order to match this approach with the one used in computing the local thermal expansion, which uses a continuous coarse graining function, we impose that the second moment of the two coarse graining functions are matched. This is achieved for r c ≈ 2σ. In the following, we use r c = 10Å, comparable to the coarse graining size σ = 5Å used to obtain the thermal internal stress in the previous section.
The probability distribution function of the yield stress, shown in Fig.3, is well fitted by an hyperbolic probability distribution function of the form: Again we use this distribution for convenience, without any particular theoretical justification. Knowing the probability distribution of the yield stress and the one of the thermal stress, and assuming that these two quantities are statistically independent, we can estimate the probability that the thermal stress triggers a local plastic event within the coarse grain size, p(σ th > σ y ): Θ(x) is Heaviside function, p(σ th , σ y ) is the joint probability of the thermal stress and yield stress within coarse grain size. We assume that the thermal stress and yield stress are   Fig.3, η and γ are fitted in Fig.2, where ∆T = 100 K independent variables, so that p(σ th , σ y ) = p L (σ th )p H (σ y ), where p L (σ th ) and p H (σ th ) are the probability distribution functions of the thermal and yield stress fitted by the expressions given above. Then the probability that a shear transformation is activated due to thermal dilation effects can be calculated as Using the parameters extracted from simulation data, we can then estimate this probability numerically. The results are shown in table I, and discussed in the following section.

V. DISCUSSION AND CONCLUSION
The probability obtained in table I is small. On the scale of a simulation box that contains a few thousands atoms (or equivalently, a few hundred possible shear transformations), the chance of actually observing a plastic event triggered by thermal stresses will be negligible.
On the other hand, this probability will be significant on the macroscopic scale of an experiment. More precisely, if we consider a sample divided in n potential shear transformation sites, the probability of observing at least one shear transformation in this sample during one temperature cycle is P act = 1−(1−p(σ th > σ y )) n , which we will refer to as the activation probability. This probability becomes larger than 90% for n ≈ 5000. This corresponds to a system of lateral size that can be roughly estimated to 10 nm, assuming that the system is made of non overlapping zones that can undergo independent shear transformations.
An important question concerns the robustness of our results with respect to the particular choice of coarse graining size. While our initial choice was motivated by the typical scale for different coarse grain size, the dash line is guided for the eye, and the value is the mean of three coarse grain size which equals 1.31 × 10 −6Å 3 . The sample quench rate is 10 9 K/s. of shear transformations in metallic glasses, it remains somewhat arbitrary. Moreover, it is seen in figure 1 (b), that the variance of the local thermal expansion depends on coarse grain size σ as a power law, so actually there is no characteristic coarse grain size for local thermal expansion. We therefore determined P act for different coarse grain size, with the idea that this activation probability is a physical value in a given sample, and should not be sensitive to the arbitrary coarse graining size σ. Since p(σ th > σ y ) 1 we can in general simplify the expression of the activation probability, P act ∼ np(σ th > σ y ), where n is the number of non overlapping coarse grain zones . For a given system, we have n ∼ 1 σ 3 . Therefore, the assumption that P act is a physical quantity implies that p(σ th > σ y ) should increase with coarse grain size, while P act ∼ p(σ th >σy) σ 3 should remain constant. In Fig.4, we show the results obtained for three different coarse graining sizes to test this hypothesis. As the coarse graining size increases, the distributions of the local yield stress and of the thermal internal stress are modified, with typically higher probabilities to obtain low values of the stress. The balance between the two evolutions leads to a nontrivial change in p(σ th > σ y ), which increases with coarse grain size. In Fig.4(d), we see that this increase leads indeed to a value of p(σ th >σy) σ 3 that just fluctuates around its average value, therefore establishing the physical character of P act . Two factors will influence the efficiency of the cryogenic rejuvenation effect. One is the number of cycles, which will progressively increase the density of events until a significant rejuvenation is achieved, bringing the system into a higher energy state. Indeed, in the experiments [5], it was observed that the sample should be trained for several thermal cycles to obtain a notable rejuvenation effect. Eventually this rejuvenation will be counterbalanced by the aging relaxation from the higher energy states, until a stationary situation is reached.
The second important factor is of course the amplitude ∆T of temperature cycles. Following equation 7, the internal thermal stress increases in proportion of the amplitude of the temperature cycle and the local yield stress is insensitive with temperature at low temperature, thus p(σ th > σ y ) would increase with the temperature amplitude, so that one could expect to amplify the effect by increasing ∆T . However, when the temperature increases, thermal relaxation will come into the picture and eventually dominate the structural evolution. Aging effects would then hide the rejuvenation effect caused by thermal expansion heterogeneity, and large thermal cycling would induce an aging effect [24,25]. Depending on specificities of each sample (thermal history, time scale of thermal relaxation) different quantitative effects may be observed. In conclusion, our calculation provides the first quantitative evidence that thermally induced internal stresses have the potential for activating mechanically the plastic activity in a glassy sample, and it provides an atomic mechanism for thermal cycling rejuvenation.