Atomistic Modelling of Thermal-Cycling Rejuvenation in Metallic Glasses

Cycling of a metallic glass between ambient and cryogenic temperatures can induce higher-energy states characteristic of glass formation on faster cooling. This rejuvenation, unexpected because it occurs at small macroscopic strains and well below the temperatures of thermally induced structural change, is important, for example, in improving plasticity. Molecular-dynamics simulations elucidate the mechanisms by which thermal cycling can induce relaxation (reaching lower energy) as well as rejuvenation. Thermal cycling, over tens of cycles, drives local atomic rearrangements progressively erasing the initial glass structure. This arises mainly from the heating stage in each thermal cycle, linked to the intrinsic structural heterogeneity in metallic glasses. Although, in particular, the timescales in MD simulations are shorter than in physical experiments, the present simulations reproduce many physically observed effects, suggesting that they may be useful in optimizing thermal cycling for tuning the properties of metallic glasses and glasses in general.

In MGs, the atomic nearest-neighbour coordination shells are diverse, and it is increasingly recognized that there is significant structural and dynamic heterogeneity over a range of length-scales. Accordingly, the strains associated with thermal expansion/contraction must have a significant non-affine component 35,36 that must be associated with the development of internal stresses. While CTC was initially aimed at rejuvenation of MGs, other work 27,37 has shown that this treatment can also induce relaxation and oscillatory (rejuvenation-relaxation) behaviour.
Relatively easy rearrangement, associated with non-directional metallic bonding, may make the effects of thermal cycling (TC) particularly significant in MGs, but the non-affine nature of thermal strains applies across all classes of non-crystalline material, and is relevant (though scarcely studied so far) for the interpretation of the evolution of their structure and properties whenever temperature changes are involved.
For MGs, molecular dynamics (MD) simulations have been useful in understanding the links between structure and dynamics, and specifically the role of structural heterogeneity 38 . MD simulations have been used to study the effects of TC on a model Lennard-Jones (L-J) glass [39][40][41] ; in this series of studies, the temperature variation with time had a triangular waveform, i.e. there was no hold-time at the upper and lower temperatures. The TC had an amplitude as high as 2-92% of the glass-transition temperature Tg. Predominantly, the effect of TC was to lower the potential energy of the glass, and the extent of this relaxation was maximum for intermediate values of the TC amplitude. In the glass formed at the lowest cooling rate and subjected to the maximum amplitude of TC, it was just possible to detect rejuvenation. The MD simulations showed that the extent of the stress overshoot on a tensile stress-strain curve was clearly greater for lower energy states. The simulations also allowed mapping of the distribution of clusters within which the non-affine displacements induced by TC are concentrated.
In the present work, we simulate TC of an actual glass, Cu50Zr50 (composition in at.%), and apply heating and cooling, at 10 14 K s -1 , that is essentially instantaneous. We find clear rejuvenation as well as relaxation, and we map the boundary between these effects. MD simulations (generally, and in the present work) involve timescales much shorter than those typical in physical experiments. Nevertheless, we show that MD can assist in understanding the link between TC conditions and property changes. In particular, the simulations can identify the relevant length-scale of heterogeneity (a question not addressable in the original work on CTC 9 ), and can identify the stage in TC mainly responsible for rejuvenation.

Results
Potential energy variation during thermal cycling. Cu50Zr50 (at.%) metallic glasses are simulated using classical MD (see Methods). Glasses formed by quenching from the liquid state at different rates Q are subjected to TC between Ta (400 K) and Tb (mostly 1 K) with hold times of ta and tb at these temperatures. The potential energy E is monitored, and Fig. 1a shows a case where its value (initially E0 at the end of the quench) is higher after each cycle. The effect of the overall processing (N cycles) is characterized by the value at the end of the holding time ta at the final temperature (400 K).
A higher Q gives an initial glass with a higher E0 (Fig. 1b). Values of E, after up to 60 cycles, are shown for ta = 10 ps and ta = 200 ps; for a given ta, independent of E0, all the potential energies converge to a common value after about 40 cycles (Fig. 1b). Following the general interpretation of the effects of energetic processing 1 , this value Es is taken to represent a steady state in which structural damage and repair rates are in balance. Comparing ta = 10 ps and 200 ps, the longer time allows for more relaxation and therefore gives a lower Es. For 10 ps, all the glasses show rejuvenation on cycling; for 200 ps, the glass formed at 10 12 K s -1 stays at roughly the same energy, and the most rapidly quenched glass (Q5 = 10 13 K s -1 ) undergoes relaxation.
For a given treatment, a more relaxed glass shows slower rejuvenation (Fig. 1b, and Supplementary Fig. 1a where the energy changes are compared from a common starting point). This is consistent with the suggestion 35 that the origin of the rejuvenation is the heterogeneity of the glass structure; a more relaxed initial glass (formed on slower cooling) would have a more uniform structure. A more direct test of the effect of uniformity is obtained by MD simulations of B2 CuZr crystal. In this case, crystallographic symmetry dictates that the thermal strains are affine.
The potential energy indeed shows no changes as a result of TC (see Fig. 6c and Fig. 7b later).
For a given initial glass, the effect of ta on the ES induced by TC is shown further in Supplementary Fig. 1b: Es depends strongly on ta and intersects the E0 values of the as-quenched glasses (Fig. 2). A given as-quenched state can thus be rejuvenated or relaxed by TC (Fig. 2), just as seen in experiment 27 . Other TC parameters being equal, there is a critical value tc (arrowed): for ta < tc there is rejuvenation, and for ta > tc relaxation. The ultimate energy change on TC (∆E s = E s − E 0 ) can be mapped as a function of ta and of the initial E0 value of a range of glassy samples ( Supplementary Fig. 2). Extrapolating the trends on the map, for a typical bulk metallic glass with cooling rate 10 4 K s -1 , the holding time at the upper temperature should be shorter than 100 s to achieve rejuvenation.
Upon TC, the model L-J glass cited earlier [39][40][41] predominantly showed relaxation. The higher of the quench rates used to form the initial glasses in the present work (>10 11 K s -1 ) overlap the quench rates applied to the L-J system, yet our glasses formed at those higher rates still predominantly show rejuvenation upon TC (Fig. 1b). This contrasting behaviour may reflect the heating and cooling rates (Qa and Qb) used in our TC; these are some two orders of magnitude higher than the TC rates applied to the L-J glass, and the consequences are explored further below.
In some of the TC simulations for the L-J glass, the upper temperature was as high as 92% of Tg, and this would also favour overall relaxation.
Structural and mechanical properties of cycled samples. Volume-average properties are compared for a glass with a given Q in its as-quenched state and after TC with two different ta values to achieve steady states that represent relaxation (ageing), and rejuvenation with respect to the as-quenched glass (these three states are indicated on Supplementary Fig. 2). In each sample, the shear and bulk elastic moduli show a distribution of values. As expected, the distributions are shifted upwards in the relaxed sample and downwards in the rejuvenated sample (Fig. 3a). Atomic pair distribution functions show minor differences between the three samples. The changes induced by rejuvenation are opposite to those induced by relaxation (Fig. 3b). The different configurations of nearest-neighbour clusters can be described as Voronoi polyhedra. For the glasses in the present work, the relative populations of the most favoured polyhedra are shown in Fig. 3c. These populations all increase on relaxation and decrease on rejuvenation. In all these respects, the relaxation and rejuvenation achieved by TC are qualitatively the same as those achieved by annealing or by forming the glass at different quench rates.
Of most practical interest is the form of the stress-strain curve for shear deformation (Fig. 3d); in a rejuvenated sample, the stress overshoot can be completely eliminated. The elimination of the overshoot is not seen for the simulated L-J glass [39][40][41] , even when formed by quenching much faster than the maximum rate in the present work. By markedly reducing, or even eliminating, the stress overshoot at yielding, rejuvenation offers the prospect of avoiding strain localization and consequent catastrophic mechanical failure, as has been explored experimentally for rejuvenation achieved by means other than TC 8,9 .
Vibrational behaviour during TC. MD simulations are of particular interest to characterize local structures and properties in MGs; studies show that there are low-frequency (soft) modes associated with a distribution of soft spots in the glass 42,43 . As described in Methods, we focus on the vibrational (phonon) density of states (VDOS) g(E) and the atomic participation ratio Pi in low-frequency (soft) modes ( Supplementary Fig. 3).
A glass subjected to TC shows an increase in ∆E = E N − E 0 , accelerating as the heterogeneity increases, before levelling off at steady state (inset, Fig. 4a). The VDOS (Fig. 4a) shows the boson peak characteristic of glasses. This peak shifts to lower energy and is more intense as the number of cycles N increases. These effects are similar to those found for heavy plastic deformation of an MG in an experimental and MD simulation study 44 , in which they were associated with rejuvenated regions at shear bands.
The atomic participation ratio Pi is a suitable basis for considering temporal and spatial correlations in the MG. The Pearson correlation coefficient p0,N of Pi, between the initial state and after N cycles at any given location shows a clear decay towards zero as N increases (Fig. 4b). The memory of the initial state is progressively lost as the sample reaches the steady state, implying that TC can induce comprehensive structural change (α relaxation) even though the temperature is always far below Tg. As the population density of soft spots (and the boson peak intensity) increase with the number of cycles, more rearrangement events are triggered during TC; thus the sample structure changes more and more rapidly. The decrease in correlation of Pi between two adjacent cycles as TC progresses (inset, Fig. 4b) reveals the accelerating increase in the heterogeneity of the MG.
The spatial distributions of Pi for the initial glass (Q1= 10 9 K s -1 ) and after the first cycle (Figs. Relaxation can also be restrained by lowering Tb, shortening tb, and using faster cooling and heating rates. With relaxation reduced by changing these parameters, more net rejuvenation is accumulated in each cycle ( Supplementary Fig. 6). Strikingly, even though the structure is already frozen at Tb, there is still some local relaxation, indicated particularly by the local kinetic energy or temperature (Fig. 5a). Restricting relaxation during the low-temperature hold, by decreasing tb, the distribution of local temperature Tb,loc is more spatially heterogeneous (Figs. 5b-d) at the beginning of the heating process; this heterogeneity could be amplified and could enhance the rejuvenation during heating. The heterogeneity of local temperature Tb,loc does not continue to decrease as tb is prolonged; finally the decrease saturates when the heterogeneity of Tb,loc is consistent with the structural heterogeneity. By introducing the quasi-equilibrium process (see Methods), the infinite-tb state can be obtained. At t b = ∞ , the degree of rejuvenation is independent of cooling rate Qb and lower temperature Tb, and correlated only with the heating process.
In a practical implementation of TC, there are likely to be hold times (ta and tb) at the upper and lower temperatures. We have shown that structural relaxation during these hold times is critical in understanding the overall effect of TC. The effects of such relaxation cannot, of course, be studied if the TC is simulated only with a triangular-wave thermal history as in Refs. [39][40][41].
The causes of rejuvenation in TC. Fixing the conditions (and therefore the likely degree of relaxation) during the upper-and lower-temperature holds, we focus on the effects of the rates of cooling and heating, Qb and Qa. The potential-energy evolutions in one cycle with slow cooling and heating (10 10 K s -1 ), and with ultra-fast cooling and heating (10 14 K s -1 ), are shown in Fig. 6a.
For low Qb and Qa, the potential-energy evolution is almost completely symmetrical on cooling and heating, and the value of E changes only when the temperature is changing. However, for high Qb and Qa, the symmetry is broken: compared to the E profile at low rate there is little change on cooling, but major changes on heating. At the start of heating at high Qa, there is a lag before E starts to rise rapidly, overshoots its initial value E0, passes through a maximum, and ends the cycle higher than E0. We note that, if the simulation were continued to longer time, E might well decay back to E0; in other words, this rejuvenation may be transient.
The evolution of E at high Qa implies that the reheating stage plays the key role in TC-induced structural change. It is of interest to examine the correlation between the changes of potential energy ∆E p = E − E 0 and kinetic energy ∆E k = E k − E k,T b . The increase in potential energy of the glass clearly lags behind the increase in kinetic energy during ultra-fast heating (Fig. 6b, red data points) while there is no similar effect on slow heating (Fig. 6b, blue data points). This suggests that, at least in these MD simulations, the glassy system during ultra-fast heating falls out of "equilibrium" by violating the equipartition theorem (ET); there must be local decoupling of the potential and kinetic energies. In the simulation, the degree of adherence to the ET can be adjusted by choosing different thermostats, and thereby different degrees of rejuvenation are achieved (Supplementary Fig. 7 and the associated discussion in Supplementary Material). By using the Langevin thermostat to force adherence to the ET, rejuvenation of the glass is suppressed (Fig. 6b, orange data points). Using the usual Nosé-Hoover thermostat, and applying TC to the B2 CuZr crystal there is again no rejuvenation (Fig. 6b, green data points). The contrasting effects of TC (rejuvenation or not) are shown in Fig. 6c. Thus, even at high Qa, rejuvenation is achieved only when: (i) there is violation of the ET under the Nosé-Hoover thermostat; and (ii) there is structural heterogeneity, as in the glass but not the crystal. This is direct evidence that an external cause, the violation of the ET during ultra-fast heating, takes on its role through an internal cause, the intrinsic structural heterogeneity, in achieving the rejuvenation in MGs.
This point can be explored further through the atomic rearrangements induced by ultra-fast heating. The local intensity of these rearrangements in a voxel i at temperature Ti is characterized by the activated rate, defined as R i ≡ −∆F i k B T i (see Methods), where ∆F i is the activation energy barrier of the voxel (proportional to the local instantaneous shear modulus 45 ), and kB is the Boltzmann constant. The barrier ∆F i is directly determined by the local structural packing, which in a MG is distributed inhomogeneously. The Ti is controlled by the fast out-of-equilibrium heating. Spatial distributions of Ri at the start of the hold at Ta after ultra-fast heating are shown for different cases in Fig. 7. For the B2 CuZr crystal, Ri is controlled by the periodic structure (Fig. 7b, corresponding to the green data points in Fig. 6b), even though there can still be considerable local temperature heterogeneity generated by the NH thermostat. For the glassy state, the Ri is still homogeneous without local temperature heterogeneity, when the Langevin thermostat is used ( Fig. 7c, corresponding to the orange data points in Fig. 6b). Thus, consistent with what has already been noted, the intense Ri that leads to rejuvenation upon reheating seems to require both structural and kinetic heterogeneity (Fig. 7d).
The rejuvenated glass (i.e. after TC using the NH thermostat) has the broadest distribution of Ri and the highest frequency of local atomic rearrangements during the hold at Ta (Fig. 8a). We calculated the local non-affinity D i to characterize the evolution of local atomic rearrangements with holding time ta. The relatively strong correlation between the local non-affinity and the local activated rate for all atoms (Fig. 8b) in the rejuvenated case emphasizes that the atomic response is mainly controlled by the local coupling between atomic structure and temperature, which is not found in the crystal state or in the glassy state with the Langevin thermostat. The spatial evolution of the non-affinity (Fig. 9) shows that, even after a long hold at Ta, the atomic rearrangements are still correlated with the local activated rate at the start of the upper-temperature hold (circled region in Fig. 7a). Over a long enough time, however, thermal relaxation must erase such rejuvenation effects, and ultimately lead to ageing.

Discussion
Thermal cycling is known to change the properties of polycrystals of metals that have low crystallographic symmetry and therefore an anisotropic single-crystal thermal expansion coefficient 46,47 . For such polycrystals, the relevant microstructural length-scale is clearly the grain diameter. It is now accepted that, well below the temperatures at which relaxation annealing would be conducted, thermal cycling, specifically CTC, is a possible method of changing the structure and properties of MGs.
From the earliest study of the effects of CTC on MGs 9, 35 , it has been supposed that heterogeneity in the glass must play a role, but, in contrast to the polycrystalline case, the relevant length-scale has not been clear. While there are now many measurements of CTC-induced property changes in MGs, these do not give any indication of the length-scale of the relevant heterogeneity. The MD simulations in the present work do permit heterogeneity to be visualized in the spatial distributions of the local values of: the atomic participation ratio Pi (Fig. 4c,d and Supplementary Fig. 4), the temperature Tb,i (Fig. 5b-d), the activated rate Ri (Fig. 7d), and the non-affinity Di (Fig. 9a-d). As already noted, these various parameters are strongly correlated.
Thus, as seen in these figures, they share a common length-scale: this is found to be 0.5-2.0 nm.
We suggest that this length-scale is relevant for understanding TC effects in MGs.
In the TC simulations for the L-J glass, the distributions of atoms showing large non-affine displacements were mapped 39,41 . In early cycles of the L-J glass formed by quenching at the highest rate, large non-affine displacements were dispersed throughout the sample, but in all other cases, the atoms showing such displacements were dispersed as small clusters. The distribution and size of these clusters are similar to, for example, those of the local participation ratio Pi (Fig.   4c,d).
This length-scale is below the resolution limit of common imaging techniques. For electron microscopy, for example, we note that the thickness of the electron-transparent thin foil is 1 to 2 orders of magnitude greater than the heterogeneity size -accordingly, the detection of structural heterogeneity is hindered by projection effects. Nanoprobe studies of MG surfaces do suggest heterogeneity on the scale down to about 1.5 nm [48], broadly consistent with the present work.
Length-scales of tens of nm 49 , or even tens of μm 18 , have also been suggested to describe the heterogeneity in MGs, but such scales cannot be studied with the simulation methods in the present work. If the length-scale of the heterogeneity giving TC effects is 0.5-2.0 nm, that is favourable for the potential application of TC to treat MG components in MEMS and NEMS devices.
An experimental study of a Cu46Zr46Al7Gd1 BMG 50 found that TC has a strong rejuvenation effect on some properties (notably the initial yield load in nanoindentation, Fy) and only a weak effect on others (e.g. hardness, H). The relaxation effect of annealing can be reversed by TC for Fy, but not for H. These effects could be consistently interpreted by considering the MG to contain soft spots in a relatively rigid matrix: TC affects mainly the soft spots, increasing both their population density and the ease of shear within them. It was inferred that a property such as Fy is strongly affected by the distribution of soft spots, while H depends mainly on the rigid matrix, which is barely affected by TC. The soft spots were associated with shear-transformation zones, and therefore directly connected with improved plasticity 50 .
The present work is remarkably consistent with this picture. As noted above, the simulations give spatial distributions of Pi, Ri and Di that all show soft spots. It is clear that TC does specifically affect the nature and population of those soft spots. Furthermore, TC is shown (Fig.   3d) to affect the onset of plastic flow much more than its continuation, consistent with the different effects measured 50 for Fy and H.
It is evident both from physical measurements 27,37 and from the simulations in the present work (Fig. 2, Supplementary Fig. 2) that the effect of TC on MGs is delicately balanced between relaxation and rejuvenation. The dynamic nature of that balance can be better understood by noting that even when, for example, there is a strong overall rejuvenation, soft spots are disappearing as well as appearing (Fig. 4c,d).
In the present work, our aim has been to advance understanding of the effects of TC on the structure and properties of MGs, and in particular how a given MG can be rejuvenated to states of higher energy or aged to lower energy. We have conducted molecular-dynamics (MD) atomistic simulations of cycling Cu50Zr50 (at.%) glasses between Ta = 400 K and Tb = 1 K. The simulations do allow exploration of the effects of a range of thermal-cycling parameters. Depending, for example, on the quench rate used to form the glass, and on the hold time at the upper temperature during cycling, the MG can undergo either rejuvenation or ageing as a result of TC.
Focusing on rejuvenation, the effects of TC are clear in progressively increased potential energy and reduced structural correlation with the as-cast glass, while the heterogeneity in the rejuvenated glass increases. As the heterogeneity increases, the rejuvenation (quantified as the potential energy) accelerates, until levelling off at a steady-state value. Rejuvenated glasses show reduced populations of the most common atom-centred clusters, reduced elastic moduli, and an increased boson peak. The effects in the simulations are in broad qualitative agreement with published experimental findings.
The simulations show that, during one cycle, the excitation (externally driven increase in 11 potential energy) occurs in large measure only on heating from the lower temperature Tb to the upper temperature Ta. We conclude that rejuvenation in TC has two causes. The first, heterogeneity in local coefficient of thermal expansion, was suggested in early work 9,35 and computed by some of us 36 , but is negligible in the present MD simulations. The second, hysteresis arising from a violation of the equipartition theorem, is found at the exceptionally high heating rates in the present simulations, but itself relies on the structural heterogeneity intrinsic to the glassy state.
The simulated Cu-Zr glasses show conclusively that cryogenic thermal cycling is capable of inducing local atomic displacements that, over progressive cycles, destroy all memory of the as-quenched structure of the glass. The properties of the TC-treated glass reach a steady state that may be relaxed or rejuvenated with respect to the as-quenched state, depending on the energy of the initial glass and the TC conditions. Simulations of the kind demonstrated in the present work, will allow exploration of the effects of TC parameters and thereby help in optimizing TC treatments.

Simulation procedures and samples. Molecular-dynamics (MD) simulations were performed
with the LAMMPS package 51 for Cu50Zr50 (at.%) metallic glass with realistic embedded-atomic-method potentials 52 . The MD time step was set to 2 fs, and a constant pressure and temperature (NPT) ensemble was employed in which the temperature and pressure were controlled by the Nosé-Hoover 53, 54 and Parrinello-Rahman 55 methods respectively. To check the thermal-expansion effect in MD, we also used a constant volume and temperature (NVT) ensemble ( Supplementary Fig. 8); this gave similar results. For temperature control in the heating step in TC, we also used the Berendsen thermostat (B) 56 , the rescaled atomic velocity approach (R), and Langevin thermostat (L) 57 . The model system contained 10,000 atoms in a cubic box with periodic boundary conditions; we compared results with a larger system, but found no size effect. We first kept the system at 2500 K (far above the melting point) for 10 ns, and then the melt was quenched at selected rates from 2500 K to Ta to prepare the initial glassy samples with different thermal histories.
Thermal cycling (TC). The different rates of quenching, Qi = 10 9 , 10 10 , 10 11 , 10 12 or 10 13 K s -1 , give virtual glassy samples in different initial states. We set the upper temperature Ta = 400 K, which is 54% of the glass-transition temperature Tg for this model system, and therefore below the temperatures normally associated with  and  thermal relaxations. Subsequent thermal cycling ( Fig. 1a) consisted of cooling to the lower temperature Tb (set at 1 K), holding at Tb, and heating back to Ta, at which all properties are measured. The cooling and heating were near-instantaneous, at~10 14 K s -1 (chosen to be faster than the quenching rates used to form the glasses). The hold times at the lower and upper temperatures are tb, set at 100 ps, and ta, treated as variable to examine possible relaxation effects. The potential energy after each cycle is monitored as a function of the number of cycles N and of ta. TC induces changes in potential energy, ∆E = E N t a − E 0 , where E 0 is the value for the as-quenched sample: ∆E > 0 represents rejuvenation (Fig. 1a) and ∆E < 0 represents relaxation (ageing). A schematic of simulation parameters in one cycle is shown in Supplementary Fig. 5.

Mechanical properties.
To measure mechanical properties, simple shear deformation with shear rate 10 7 s −1 was employed to obtain the strain-stress curves (Fig. 3d) of each type of sample at 400 K. Ten independent loading runs were employed for each sample to ensure statistical validity.
Three types of samples ( Supplementary Fig. 2, the three white points in the map): as-quenched (Q5), aged (ta = 200 ps), and rejuvenated (ta = 10 ps) were selected to investigate the influence of TC.
Quasi-equilibrium state ( t b = ∞ ). As in a crystalline solid, in the glass at low temperature structural relaxation is suppressed. The thermal motion in the glass can then be regarded as atomic vibration around the inherent structure 58 . At the low temperature in our case, Tb = 1 K, this atomic vibration can be approximated as a harmonic oscillator. In a classical system, the equipartition theorem 59 should work well at equilibrium, and thus for the quasi-equilibrium glassy state, we have where Ek is the kinetic energy of the system, and ∆E p is the difference between the potential energy at temperature Tb and the potential energy of the inherent structure.
In the MD simulations, first we used a conjugate-gradient algorithm to obtain the inherent structure. We injected kinetic energy Ek = 2kB into that structure, and used the microcanonical ensemble (NVE) to relax the system. Accordingly, the system is in the quasi-equilibrium state when E k = ∆E p = 1k B . We have checked whether the sample starts from the quasi-equilibrium state in each cycle. In the simulation, the states obtained by TC are insensitive to the quench rate Qb, holding time tb and minimization method. For the NVE relaxation, we used 2 ps in the simulation, and longer relaxation times do not change the results.
Local property characterization. The local parameters centred in 30×30×30 voxels (cubic, of side length~1 nm) were calculated by a coarse-grained method.

Vibrational density of states (VDOS) and atomic participation ratio (Pi). The dynamic matrix
(DM) was calculated from energy-minimized samples, and the DM was diagonalized to obtain the vibrational density of states and the vibration mode. To extract the contribution of the Debye vibration mode, we used the reduced VDOS, VDOS/ω 2 . The boson peak in the low-frequency regime is shown in Fig. 4a. The participation ratio (PR) was calculated following the method used in Ref. [43], using a low 1% vibration mode. The PR of atom i is given by where e ω i is the vector polarization at atom i at vibration frequency ω. The local PR of voxel i (P i ) is given by P i = 1 N i j∈N i P j , where N i is the number of atoms in voxel i, and Pj is the PR of atom j. In glasses, the low-frequency vibration modes are localized, and the related atomic participation ratio Pi can be regarded as an indicator of local structural instability. Thus, the evolution of Pi could provide deep understanding about the local structural responses during rejuvenation.
Local activated rate. We define the local activated rate Ri of voxel i as where ∆F i is the local energy barrier, and T i is the local temperature, of voxel i.
In the 'shoving model' for molecular rearrangements in liquids 45 , the relaxation energy barrier is dominated by elastic energy, and the relaxation between energy barrier and elastic energy is even true at local scale 60 . Following the formula by Wang et al. 61 , here we assume that the formula still works at the local scale. Then the local energy barrier ∆F i can be expressed as: where G ∞ i and B ∞ i are the local instantaneous shear and bulk moduli, respectively; Vi is the coarse-grained volume of, and Ni is the number of atoms contained in, voxel i. Following the method used in Ref. [62], we use the local Cauchy-Born modulus to approximate the instantaneous modulus. The local temperature Ti is defined by: where E k i is the local kinetic energy, v j is the velocity of atom j and m j is the mass of atom j.
Local non-affinity. Following the work of Falk & Langer 63 , the non-affinity of atom j can be defined as , and r jk = r k − r j is the position vector of atom j at time t. We take Dj, the square root of D j 2 , as the non-affinity of atom j, and local non-affinity of voxel i is given by We take the beginning of the holding time at Ta as the reference time (ta = 0), as shown by the circled region in Fig. 7a, and investigate the evolution of non-affinity with ta.

Pearson correlation coefficient.
This coefficient is employed to characterize the correlation between X and Y as where E[ − ] is the expectation.     The corresponding 3D spatial distributions of Tb,loc at the end of the low-temperature hold.        reheating/cooling rate Qa = Qb = 2.5×10 14 K s -1 , tb = 200 ps, ta = 4 ps, Tb = 1 K and Ta = 400 K.

The evolution of soft spots in TC
In TC, the soft spots in metallic glasses can be activated, and during the atomic rearrangement they can be produced or annihilated. There are examples of both disappearance and generation of soft spots (dashed circles in Fig. S4), but overall an increase in their population. With increasing number of cycles the increased population reaches saturation. For the initial few cycles (Fig. S4), most soft spots are still stable, indicating that the rejuvenation is quite localized. In contrast, when the population density of soft spots reaches saturation, the spatial distribution is noticeably shuffled; this is caused by a sequence of interactions between soft spots. ps, ta = 4 ps, Ta = 400 K and Tb = 1 K. 41 As shown in Fig. S7a, the degree of rejuvenation is highly dependent on the choice of thermostat.
In MD simulations, the thermostat is intended to mimic the processes of 'equilibration' in the physical world. It is difficult to say which of various possible thermostats best represents the physical world in a highly non-equilibrium heating process. However, the basic principle should be consistency between MD and the physical world with regard to different degrees of rejuvenation.
When the distribution of local kinetic energy at the start of the first heating stage is more heterogeneous (Fig. S7b, d), the degree of rejuvenation during heating is greater. For the Nosé-Hoover thermostat, the heterogeneity is greatly enlarged after heating while, in contrast, the Berendsen and the rescaled thermostats maintain the initial heterogeneity. Since the Langevin thermostat controls the temperature by a stochastic equation, the heterogeneity of local kinetic energy is maintained at the lowest level, and there is no correlation with the initial state. The mean value of the local temperature is identical in all cases. Thus, there is a special role for the heterogeneity of local kinetic energy, which, in addition to the glassy structure itself, can also drive the rejuvenation or relaxation.

The effect of bulk thermal expansion in simulations
Due to the heterogeneous structure of metallic glasses, non-uniform thermal expansion has been considered as the origin of rejuvenation in cryogenic thermal cycling (CTC). In our MD simulations and previous work [11], we find that the probability of activation by the local thermal expansion mismatch is quite small. If we suppress bulk thermal expansion during TC by fixing the sample volume, the efficiency of rejuvenation does not change significantly (Fig. S8), that is to say, bulk thermal expansion is unnecessary for CT induced rejuvenation. There is a role for structural heterogeneity, but not for overall volume change.