# Understanding the mechanisms of amorphous creep through molecular simulation

See allHide authors and affiliations

Edited by Sharon C. Glotzer, University of Michigan, Ann Arbor, MI, and approved October 27, 2017 (received for review May 24, 2017)

## Significance

The individual and collective molecular displacements in an amorphous solid undergoing plastic deformation are simulated by an atomistic method that allows incremental motions to be observed over a time window of fractions of seconds. Because the timescale matches well with the conditions of experimental measurements, simulation details provide dynamical evidence for the fundamental mechanisms of amorphous creep. In particular, knowledge of the interplay between diffusion (flow) and mechanical deformation processes enables us to explain the stress and temperature behavior of the experimental data as well as the validity of model descriptions of molecular mechanisms in terms of spatially and temporally heterogeneous fluctuations.

## Abstract

Molecular processes of creep in metallic glass thin films are simulated at experimental timescales using a metadynamics-based atomistic method. Space–time evolutions of the atomic strains and nonaffine atom displacements are analyzed to reveal details of the atomic-level deformation and flow processes of amorphous creep in response to stress and thermal activations. From the simulation results, resolved spatially on the nanoscale and temporally over time increments of fractions of a second, we derive a mechanistic explanation of the well-known variation of creep rate with stress. We also construct a deformation map delineating the predominant regimes of diffusional creep at low stress and high temperature and deformational creep at high stress. Our findings validate the relevance of two original models of the mechanisms of amorphous plasticity: one focusing on atomic diffusion via free volume and the other focusing on stress-induced shear deformation. These processes are found to be nonlinearly coupled through dynamically heterogeneous fluctuations that characterize the slow dynamics of systems out of equilibrium.

Deformation and flow are fundamental in the rheological behavior of many materials (1⇓⇓–4). For the molecular understanding of plastic response under stress (creep), a standing challenge is to know the details by which the constituent atoms rearrange themselves individually and collectively in a local environment of stress and temperature. Creep experiments have been extensively reported on amorphous materials, including metallic glasses and colloidal systems (5⇓⇓⇓–9). Yet, the relationship between stress and temperature effects on creep and the underlying microscopic processes remains an open question. Theoretical models have been proposed to describe the mechanisms of molecular deformation and flow responsible for amorphous plasticity. Spaepen (10) considered the distinction between homogeneous and inhomogeneous flows in metallic glasses and introduced the concept of the local free volume as an order parameter. In this view, local strain production and dissipation are assumed to be associated with individual atomic jumps. Argon (11) proposed a plastic deformation model of metallic glasses based on the notion of local shear transformations driven by stress and in the presence of thermal fluctuations. The atoms participating in such processes essentially undergo an inelastic shear deformation. Falk and Langer (12) later introduced the term shear transformation zone (STZ) in interpreting simulation results of viscoplastic deformation of amorphous solids. The STZ theory was further extended to be capable of describing temperature- and rate-dependent amorphous plasticity (13, 14). The term STZ has become widely adopted in studies of amorphous materials. (4, 8, 15⇓⇓⇓⇓–20).

We seek to identify the elementary processes of deformation and flow in amorphous creep through atomistic simulation. A bottleneck well-known in the literature is that the temporal scales relevant to creep are beyond the reach of traditional molecular dynamics (MD). To compensate for its inherently microscopic timescales, MD simulations of creep had to resort to extreme conditions of stress, temperature, and strain rate (21, 22). Here, we implement a metadynamics formulation that allows transition-state trajectories to be generated on appreciably longer timescales on the order of fractions of seconds. By analyzing the distributions of atomic strain and nonaffine particle displacement, we observe the effects of stress and temperature on the evolution of activated states at the microscale. We find that the processes of single-particle diffusion and of shear deformation of small clusters of particles are both active in steady-state creep; in particular, their interplay gives rise to a characteristic upturn behavior of stress effects on creep rate. The simulation results also support a mechanism map showing a regime of low stress and high temperature where diffusional creep dominates and a high-stress regime governed by shear deformation creep. Our findings suggest that, while the two original models of amorphous plasticity (10, 11) are complementary in their individual focus, their combined effects need to be analyzed using more fundamental theories capable of treating the effects of nonlinear feedback.

## Simulation Methods and Model

We consider a model Cu_{50}Zr_{50} metallic glass system in two dimensions. The atoms interact through a Lennard–Jones potential (23, 24), which has been used to study plastic deformation (25) and thermally activated flows in metallic glass (26). Two amorphous structures with different sizes, 39.5 × 19.7 nm (containing 10,000 atoms) and 62.5 × 31.2 nm (25,000 atoms), are prepared by quenching from a high-temperature liquid state (*Materials and Methods* has additional details). We apply uniaxial tensile stress to the system and follow the procedure described in *Materials and Methods* to simulate the time evolution of creep. The algorithms allow us to investigate creep strain evolution and the corresponding molecular mechanisms at various stresses and temperatures, in particular at low stress and slow creep rate. To quantify the atomic-level plastic deformations in creep, we compute two local strains: deviatoric strain *η*_{Mises}. In addition, we compute nonaffine displacement as a measure of single-particle dynamics. Detailed descriptions of atomic strains and displacement are discussed in *Materials and Methods*. We find negligible system size effects on the results; therefore, our discussions will refer to the larger system unless explicitly stated otherwise.

## Creep Curves

Fig. 1 shows the time development of the system strain determined by metadynamics simulation. The creep curve shows the classical behavior of three stages of strain evolution. The initial period of strain nucleation and distribution, consisting of a steep increase followed by a gradual approach to saturation, is known as primary or transient creep. The secondary stage of steady-state creep is a period of linear strain increase in time. The extent of this stage depends on the combination of applied stress and system temperature. For relatively low stresses, the secondary stage may have a considerable extent before the onset of tertiary creep, where the strain rate increases without apparent limit (Fig. 2). In the case of Fig. 1, the onset of structural instability at a strain level of ∼2.5% is readily observed. Notice that the relevant timescale of creep is of the order of fractions of seconds, comparable with those of laboratory measurements and beyond the capabilities of traditional MD.

Fig. 2 shows a series of creep curves spanning the stress range from 42 to 258 MPa simulated at a temperature of 0.68*Supporting Information*) (8, 27).

## Nonaffine Displacement and Deviatoric Strain Distributions (Dynamical Heterogeneities)

Fig. 2 indicates that the applied stress has a significant effect on the time evolution of the creep (system-level) strain. One can ask for the atomic rearrangements and displacements associated with the observed deformation. In Fig. 3, we show how the probability distributions of nonaffine particle displacement

The nonaffine displacement *η*_{Mises} and *η*_{Mises} is in *Supporting Information*). Fig. 4 shows the variation with stress of the distributions of *Inset* shows that the number of atoms participating in large strain deformation indeed increases sharply with increasing stress. Corresponding to the stress regime separation noted in Fig. 3, an appreciable increase in the number of atoms involved in shear transformations is seen starting at a stress of 143 MPa. This suggests that stress plays a fundamental role in the physical manifestation of deformation and flow behavior in driven systems.

To see what molecular mechanisms could be associated with Figs. 3 and 4, we show in Fig. 5 the local space–time distributions. The spatial maps are resolved at the nanoscale and temporally measured over a time interval of 0.01 s. Maps of nonaffine displacements *A* and *B*. Correspondingly, maps of *C* and *D*. Fig. 5 *A* and *C* refers to low (64 MPa) stress, and Fig. 5 *B* and *D* for high (258 MPa) stress.

In Fig. 5*A*, the features of the nonaffine displacement map can be essentially attributed to thermally activated atomic diffusion, as initially proposed by Spaepen (10). It could serve as the baseline for measuring responses through single-particle displacements at higher stress. Fig. 5*B* indicates that the displacement magnitude of the sites of high mobility is about a factor of three larger than the displacements at low stress. Moreover, we see the presence of bidirectional flow (red–blue interface in Fig. 5*B*), as several such active zones appear along the 45° direction. When looking at Fig. 5*D*, these activated events appear as local strain bursts, a striking effect of stress-induced deformation. We regard this to be a specific simulation result, obtained at a physically meaningful timescale, that illustrates the synergistic effects of stress activation.

In Fig. 5*C*, we see little or no local shear deformation activity, which suggests that thermally diffusional rearrangements at the temperature of the simulation produce relatively small local plastic strain. These weak events seem quite random in space and do not align in any particular direction. We conclude that, at low stress, creep proceeds mostly by single-particle diffusion rather than cooperative atom deformation. However, at high stress, Fig. 5 *B* and *D* implies that diffusive and deformation processes are both activated and coupled through heterogeneous stress fields in the driven system.

We apply the above findings to predict the variation of creep rate with tensile stress, a behavior of fundamental interest in the mechanics and physics of creep (30, 31). At very low stress, we expect the effects of thermal activation to dominate, which should lead to a mild increase of creep rate with stress. As stress increases, both flow and deformation mechanisms contribute smoothly until the stress exceeds a threshold, at which point thermal and stress activations accelerate. The significant increase in shear transformation deformation events, shown in Fig. 4, *Inset* and the gap noted in Figs. 3 and 4, are hallmarks of a threshold (cross-over) behavior that is widely observed in experiments (Fig. 6*B*). We consider this stress threshold to reveal a regime of nonlinear response, which has not been previously interpreted at the molecular level.

## Creep Rate Upturn

The simulated creep curves allow us to determine a steady-state creep rate at each stress. Plotting the results in Fig. 6*A*, we see a bimodal behavior in the monotonic variation of creep rate

To compare our results with experiments, the simulation data are shown again in Fig. 6*B* along with two sets of measurements (9, 33). Reduced stress

## Creep Mechanism Map

The creep simulations that we have conducted cover three temperatures, 0.57, 0.68, and 0.91

Looking at the simulation data across the temperature range, one can visualize a mechanism boundary separating the map into two domains, where either diffusion or shear deformation predominates. To motivate such a boundary, we recall well-known empirical expressions based on transition-state theory for the limiting behavior of strain rate in amorphous systems. At low stress, the creep rate *Q*_{d} and *V*_{d} are the activation energy and activation volume for diffusion, respectively. This is derived based on the free volume diffusion model. At high stress, a more appropriate form is (8, 11)**1** equal to Eq. **2** and treating these activation parameters numerically, we find the locus of (σ, *T*) values for the mechanism boundary to be essentially a straight line: *Supporting Information* has details). The red dashed curve in Fig. 7 is drawn having this functional form and adjusted to separate the groups of squares and circles. Thus, the simulation data are able to accommodate a boundary that would have the two limiting behaviors described by Eqs. **1** and **2**.

For experimental validation, we show in Fig. 7 mechanical deformation measurements on a five-component metallic glass (9). Except for an overall shift, the experimental and simulation results both display the same predominant domains of diffusional and deformational creep. One can offer only qualitative reasons why a comparison between this simulation and experiment should be considered meaningful. For example, differences between the simulation system and the experimental samples can be quite significant. These could arise from their respective preparations, such as cooling rate and casting defects and surface polish imperfections. Another factor is that an idealized 2D binary mixture model is being compared with measurements on a 3D five-component material. Although a 2D model could be capable of capturing aspects of thermally activated diffusion and stress-induced shear transformation processes, one could argue that the system dimensionality could affect activation energy barriers of elementary processes, which in turn, could influence the creep rate. In addition, the cooling rate in sample preparation will affect the underlying potential energy landscape, and therefore the creep rate response. Indeed, we have performed additional simulations at a cooling rate three orders of magnitude higher and found the mechanism boundary to shift downward in a manner that brings experiments and simulations in Fig. 7 closer together (a discussion of cooling rate is in *Supporting Information*). Other than sample preparation, the shift of the mechanism boundary also may be rationalized by noting that different spatial correlation lengths are involved in single-particle jumps vs. collective distortions of particle clusters, the former being more short-ranged and therefore less cooling-rate sensitive than the latter (35).

## Discussion

In this work, we relate the simulation results to the mechanisms proposed by Spaepen (10) and by Argon (11), because each provides a clear-cut physical description for identifying the molecular processes governing creep under varying stress and temperature conditions. An appropriate next step would be to connect the simulation results to theoretical analysis, such as the STZ concept of plastic deformation in amorphous materials. The term STZ, initially conceived to describe a transient flow defect, has evolved as a statistical thermodynamics-based theoretical description of system-level response in a thermal and stress environment. The original theory (12) has been reformulated (14), reviewed (36), and updated (37). At the this time, STZ seems to play two different roles: a mechanism combining the ideas of Spaepen (10) and Argon (8, 11, 38) or a theoretical framework for interpreting experiments or simulations (39, 40).

In discussing Fig. 6, we have regarded the creep rate upturn as signifying the synergistic actions involving single-particle diffusion and shear-induced deformation. One can imagine quantifying this interpretation by using STZ theory to explicitly calculate the creep rate response to stress. Recently, Langer (37) has reported a study of a similar problem, in which the density of STZs and an effective thermodynamic temperature were introduced as dynamical variables in a set of coupled equations of motion with a stress-dependent deformation rate.

We believe that the increase of the stress exponent from approximately two to five delineates a regime of nonlinear response in the rheological behavior of amorphous materials, a consequence of the nonlinear coupling between atomic diffusion and shear-induced deformation. More generally, it indicates an interplay between thermal and stress-activated processes, involving loading effects, thermal noise, and stress-induced fluctuations (41).

In summary, we present findings to quantify the molecular mechanisms of steady-state creep in a model metallic glass. At low stress and high temperature, the dominant mechanism is observed to be thermally activated particle flow, while at high stress, the mechanism is a more complex process of stress-induced enhanced local shear deformation and atomic diffusion. Taking these processes together leads to an interpretation of the experimentally observed stress and temperature behavior of amorphous creep as well as a unifying picture of single-particle diffusion and collective atom rearrangements. This perspective motivates revisiting the existing notion of dynamical heterogeneities (42⇓–44) and a variant of self-organization in slowly driven threshold systems (45) in the spirit of assessing various theoretical frameworks, such as the STZ theory (37), expanded mode-coupling formalism (46), mean field approaches with weakening mechanism (47), and time-dependent transition theory modeling (48).

## Materials and Methods

### Model Metallic Glass.

To prepare the metallic glass thin film configuration for atomistic simulations, we quench from a high-temperature, well-equilibrated liquid state at 0.14 K/ps at zero pressure using MD with periodic boundary conditions in all directions. We observe an inflection point in the volume–temperature curve at 440 K, which is regarded as the glass transition temperature *T*_{g} (*Supporting Information*). To model a uniaxial stress experiment, two free surfaces in the y-direction are created by removing periodic boundary conditions, so that the stress is free in that direction. We perform another 400-ps MD simulation to relax the system to zero average stress. The mechanical properties of the thin film, such as elastic modulus and yield strength, are then characterized (*Supporting Information*).

### Metadynamics Simulation of Creep.

To simulate creep using a metadynamics algorithm, we apply a prescribed uniaxial tensile stress to the system and execute the following steps.

*i*) Perform energy minimization on the relaxed system to bring it to the nearest local energy minimum.*ii*) Apply the autonomous basin climbing (ABC) algorithm (49, 50) to obtain the transition-state pathway and determine the neighboring local energy minimum state.*iii*) Compare the internal stress of the new state with the prescribed tensile stress. If the two stresses deviate by more than 1%, perform step*iv*; otherwise, go back to step*ii*.*iv*) Perform cell relaxation in the presence of the external stress. The atoms are rescaled to new positions whenever the size of simulation cell is changed, and the final configuration converges to a new local minimum.

The output of ABC is a set of transition-state pathway trajectories, each being an ordered sequence of energy minima and saddle points. The system evolution is then determined by examining the newly sampled configurations of the local energy minima. The activation time of each evolution step is estimated through transition-state theory, ^{12} s^{−1}.

### Atomic-Level Strain and Nonaffine Displacement.

We consider two atomic-level strains, the deviatoric strain *η*_{Mises}, as measures of local plastic deformation (12, 52). Imagine that a region surrounding an atom undergoes a strain deformation during a time interval *η*_{Mises} in two dimensions is computed by*η*_{Mises} is directly derived from 𝐉, while *η*_{Mises} and *n* neighboring atoms. Both can serve as an indicator for cooperative local rearrangement.

The nonaffine displacement is useful for tracking single-atom dynamics. For a time interval (

## Acknowledgments

We thank A. S. Argon, J. S. Langer, and K. Kamrin for discussions, and the Massachusetts Institute of Technology (MIT) Nuclear Science and Engineering Department for computational resources. P.C. was supported by US Department of Energy (DOE) Grant DE-NE0008450. M.P.S. was supported by National Science Foundation CAREER Grant DMR-1654548. S.Y. was supported by DOE-Basic Energy Sciences Grant DE-SC0002633 and the Kuwait-MIT Center Signature Project. S.Y. is a founding member of the MIT Concrete Sustainability Hub.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: syip{at}mit.edu.

Author contributions: P.C. and S.Y. designed research; P.C. performed research; P.C., M.P.S., and S.Y. analyzed data; and P.C., M.P.S., and S.Y. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1708618114/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- ↵
- Johnson WL,
- Lu J,
- Demetriou MD

- ↵
- ↵
- Schall P,
- Weitz DA,
- Spaepen F

- ↵
- Krisponeit JO, et al.

- ↵
- Huang Y,
- Shen J,
- Chiu Y,
- Chen J,
- Sun J

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Falk ML,
- Langer JS

- ↵
- Falk M,
- Langer J,
- Pechenik L

- ↵
- Bouchbinder E,
- Langer J

- ↵
- ↵
- Tanguy A,
- Leonforte F,
- Barrat JL

- ↵
- Pan D,
- Inoue A,
- Sakurai T,
- Chen M

- ↵
- Cao P,
- Park HS,
- Lin X

- ↵
- Greer A,
- Cheng Y,
- Ma E

- ↵
- Cao P,
- Lin X,
- Park HS

- ↵
- ↵
- Sentjabrskaja T, et al.

- ↵
- ↵
- Deng D,
- Argon A,
- Yip S

- ↵
- ↵
- ↵
- ↵
- ↵
- Weeks ER,
- Crocker JC,
- Levitt AC,
- Schofield A,
- Weitz DA

- ↵
- Cottrell AH

- ↵
- Nabarro FRN,
- De Villiers F

- ↵
- Boyle JT,
- Spence J

- ↵
- Nieh T,
- Wadsworth J

- ↵
- Klueh R

- ↵
- ↵
- ↵
- Langer J

- ↵
- Hufnagel TC,
- Schuh CA,
- Falk ML

- ↵
- Langer J

- ↵
- Kamrin K,
- Bouchbinder E

- ↵
- ↵
- ↵
- ↵
- Becker T,
- Smith JC

- ↵
- Jensen HJ

- ↵
- Gruber M,
- Abade GC,
- Puertas AM,
- Fuchs M

- ↵
- ↵
- Fan Y,
- Yildiz B,
- Yip S

- ↵
- ↵
- Cao P,
- Li M,
- Heugle RJ,
- Park HS,
- Lin X

- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Engineering