# Spin-circuit representation of spin-torque ferromagnetic resonance

Kuntal Roy<EMAIL_ADDRESS>Department of Electrical Engineering and
Computer Science
Indian Institute of Science Education and Research Bhopal, Bhopal, Madhya
Pradesh 462066, India

###### Abstract

Spin-torque ferromagnetic resonance (ST-FMR) particularly using magnetic
insulators and heavy metals possessing a giant spin Hall effect (SHE) has
gotten a lot of attention for the development of spintronic devices. To devise
complex functional devices, it is necessary to construct the equivalent spin-
circuit representations of different phenomena. Such representation is useful
to translate physical equations into circuit elements, benchmarking
experiments, and then proposing creative and efficient designs. We utilize the
superposition principle in circuit theory to separate the spin Hall
magnetoresistance and spin pumping contributions in the ST-FMR experiments. We
show that the proposed spin-circuit representation reproduces the standard
results in literature. We further consider multilayers like a spin-valve
structure with an SHE layer sandwiched by two magnetic layers and show how the
corresponding spin-circuit representation can be constructed by simply writing
a vector netlist and solved using circuit theory.

In ferromagnetic resonance (FMR) experiments, Silsbee, Janossy, and Monod
(1979); Heinrich  _et al._ (1987); Guan  _et al._ (2007); Wang, Ramaswamy, and
Yang (2018) a magnetization precesses around an effective magnetic field and
the rotation is sustained by a transverse ac field acting as negative damping.
Landau and Lifshitz (1935); Gilbert (2004) In spin-torque ferromagnetic
resonance (ST-FMR) experiments, an in-plane alternating current in spin-orbit
materials drives the magnetization precession in an adjacent magnet via direct
spin Hall effect (SHE) Dyakonov and Perel (1971); *RefWorks:764;
*RefWorks:1198; Murakami, Nagaosa, and Zhang (2003); *RefWorks:765;
*RefWorks:902; Valenzuela and Tinkham (2006); Roy (2014) and it pumps spins
into the spin-orbit layer detected as a charge voltage/current via inverse
spin Hall effect (ISHE). Wang  _et al._ (2006); Jiao and Bauer (2013); Azevedo
_et al._ (2005); Saitoh  _et al._ (2006) Spin pumping Tserkovnyak, Brataas,
and Bauer (2002a); *RefWorks:1041; Tserkovnyak  _et al._ (2005) is the
reciprocal phenomenon Brataas  _et al._ (2012) of spin-torque Slonczewski
(1996); *RefWorks:155 according to Onsager’s reciprocity. Onsager (1931a);
*RefWorks:1293 Although ST-FMR was developed first for metallic magnets, Liu
_et al._ (2011); Wang  _et al._ (2014); Zhang  _et al._ (2015) later it has
been shown that ST-FMR is also possible for magnetic insulators e.g., yttrium-
iron-garnet (YIG), Jungfleisch  _et al._ (2017) which is promising for the
development of spintronic devices. Spin Hall magnetoresistance (SMR) Chen  _et
al._ (2013); Nakayama  _et al._ (2013); Chiba, Bauer, and Takahashi (2014);
Chiba  _et al._ (2015); Schreier  _et al._ (2015); Jungfleisch  _et al._
(2015); Sklenar  _et al._ (2015) i.e., the dependence of electrical resistance
of spin-orbit material on the magnetization direction has been observed for
magnetic insulators due to the simultaneous action of SHE and ISHE, and it is
modulated by the magnetization direction via the interfacial spin transfer.
The fundamental difference of SMR with other magnetoresistances (anisotropic
magnetoresistance (AMR), Thomson (1856) giant magnetoresistance (GMR), Baibich
_et al._ (1988); *RefWorks:433 tunneling magnetoresistance (TMR) Julliere
(1975); *RefWorks:76; *RefWorks:33; *RefWorks:74) is that current does not
need to pass through the magnet. Both longitudinal SMR and spin pumping
contribute to the dc voltage/current detected in ST-FMR experiments.

Here we develop the spin-circuit representation of ST-FMR comprising both SMR
and spin pumping. Circuit theory (Kirchhoff’s current and voltage laws, KCL
and KVL, originating from the conservation of charge and energy, respectively)
has been tremendously successful for the development of the transistor-based
technology. Rabaey, Chandrakasan, and Nikoliç (2003) The physical equations
are translated first into circuit elements for simplified understandings and
it helps us to apply the principles of circuit theory e.g. superposition
principle. Then we further analyze the circuit in this organized form, solve
it programmatically, benchmark experiments, and can develop complex creative
designs for applications. Rabaey, Chandrakasan, and Nikoliç (2003) Due to the
advent of new materials and phenomena in the field of spintronics, it needs to
connect different spin-circuit representations of different phenomena and this
is necessary for the commercial developments e.g., SPICE (Simulation Program
with Integrated Circuit Emphasis). hsp  Such spintronic circuit models apply
to quantum transport Sah (1971) as well. The voltages and currents at
different nodes are of 4-components (1 for charge and 3 for spin vector) and
the conductances are $4\times 4$ matrices ($c$-$z$-$x$-$y$ basis). Some spin-
circuit representations that have been developed earlier are for ferromagnet
(FM), normal metal (NM), FM-NM interface, spin Hall effect, and spin pumping.
Brataas, Bauer, and Kelly (2006); Srinivasan  _et al._ (2014); Hong, Sayed,
and Datta (2016); Roy (2017a)

Figure 1 shows the schematic diagram of the ST-FMR setup. A longitudinal ac
current in the SHE layer having a length $l$, width $w$, and thickness $t$ in
the presence of a dc magnetic field drives the magnetization in the adjacent
magnetic layer into precession via SHE and utilizing a bias-tee we can detect
the dc voltage/current in the longitudinal direction due to SMR and spin
pumping.

Figure 1:  Schematic diagram of ST-FMR setup and axis assignments. With the
help of a bias-tee, simultaneous ac transmission and dc voltage/current
detection is possible in ST-FMR. A longitudinal ac current in $x$-direction
$J_{c}$ across the surfaces marked by 1 and 2 of the SHE layer generates a
transverse spin current density $J_{s}$ and it can drive the magnetization
$\mathbf{m}$ into precession. $H$ is the applied dc magnetic field, $H_{ac}$
is the rf driving field produced from $J_{c}$ with a phase lag $\delta$,
$\omega$ and $\theta$ are the precession frequency and angle, respectively.
Spin potentials are developed across the surfaces marked by 3 and 4. Using a
bias-tee, the dc voltage $V_{dc}$ can be detected and it comprises of both the
contributions of longitudinal SMR and spin pumping due to ISHE. Also, charge
potentials are developed across the surfaces marked by 5 and 6 due to
transverse SMR. In the axis assignments, standard spherical coordinate system
is used with azimuthal angle $\phi$ and polar angle $\theta$ for magnetization
$\mathbf{m}$, and polar angle $\theta_{H}$ for the dc field $H$.

To develop the spin-circuit representation we can utilize the Ohm’s law with
spin-orbit interaction (including both SHE and ISHE) Chen  _et al._ (2013)

$\left[{\begin{array}[]{c}\mathbf{J^{c}}\\\ \mathbf{J^{z}}\\\
\mathbf{J^{x}}\\\ \mathbf{J^{y}}\\\
\end{array}}\right]=\sigma\left[{\begin{array}[]{cccc}1&\theta_{sh}\mathbf{\hat{z}}\times&\theta_{sh}\mathbf{\hat{x}}\times&\theta_{sh}\mathbf{\hat{y}}\times\\\
\theta_{sh}\mathbf{\hat{z}}\times&1&0&0\\\
\theta_{sh}\mathbf{\hat{x}}\times&0&1&0\\\
\theta_{sh}\mathbf{\hat{y}}\times&0&0&1\\\
\end{array}}\right]\left[{\begin{array}[]{c}-\mathbf{\nabla}V^{c}\\\
-\mathbf{\nabla}V^{z}\\\ -\mathbf{\nabla}V^{x}\\\ -\mathbf{\nabla}V^{y}\\\
\end{array}}\right]$

together with the charge and spin diffusion equations $\mathbf{\nabla\cdot
J^{c}}=-\sigma\nabla^{2}V^{c}=0$ and $\mathbf{\nabla\cdot
J^{x,y,z}}=-\sigma\nabla^{2}V^{x,y,z}=-(\sigma/\lambda^{2})V^{x,y,z}$,
respectively, where $\mathbf{J^{c}}$ and $\mathbf{J^{x,y,z}}$ are current
densities, and $V^{c}$ and $V^{x,y,z}$ are spin potentials, $\sigma$ is the
conductivity, $\lambda$ is the spin diffusion length, and $\theta_{sh}$ is the
spin Hall angle of the SHE layer.

Since we are not considering any charge component in the magnetic insulators,
the generation of inverse spin Hall voltage is in the transverse direction due
to the symmetry involved in the system, and unlike charge pumping, Brouwer
(1998); *RefWorks:1324 a precessing magnet injects a _pure_ spin current into
surrounding conductors, Tserkovnyak, Brataas, and Bauer (2002a);
*RefWorks:1041; Tserkovnyak  _et al._ (2005) we consider a reduced 3-component
version of the spin-circuit and show that it can reproduce the established
expressions in literature for the dc voltage/current generated in ST-FMR
experiments comprising both SMR and spin pumping. We further employ such
formalism for multilayers like a spin-valve structure with an SHE layer
sandwiched by two magnetic layers. Chen  _et al._ (2013) We show how we can
simply write a vector netlist comprising the conductances and voltage/current
sources so that we can solve for the voltages/currents at the different nodes
_without_ invoking any boundary condition using a circuit solver. For complex
structures, derivation of any analytical expression becomes tedious and this
signifies the prowess of the spin-circuit approach.

Figure 2(a) shows the spin-circuit representation of ST-FMR that constitutes
the spin-circuit representations of SHE layer, the charge current dependent
spin current sources due to the charge to spin conversion in the SHE layer,
and spin pumping. The instantaneous spin pumping can be represented by
$\left[I_{sp}\right]=\left[G^{\uparrow\downarrow}\right]\left[V_{sp}\right]$,
where $\left[V_{sp}\right]=(\hbar/2e)\,[0\;0\;1]^{T}$,
$[G^{\uparrow\downarrow}]$ is the interfacial _bare_ complex spin mixing
conductance between the magnetic layer and the SHE layer represented in
$(\mathbf{m},\mathbf{m^{\prime}},\mathbf{m}\times\mathbf{m^{\prime}})$ basis
(where ′ denotes the time derivative) as

$\left[G^{\uparrow\downarrow}\right]=\left[\begin{array}[]{crr}0&0&0\\\
0&2G_{r}&2G_{i}\\\ 0&-2G_{i}&2G_{r}\end{array}\right],$ (1)

$\left[I_{sp}\right]=(\hbar/2e)\,[0\;2G_{i}\;2G_{r}]^{T}$ in
$(\mathbf{m},\mathbf{m^{\prime}},\mathbf{m}\times\mathbf{m^{\prime}})$ basis,
and $2G_{r(i)}$ is the real (imaginary) part of the _bare_ complex spin mixing
conductance $G^{\uparrow\downarrow}$. The _bulk_ diffusion in the SHE layer
can be represented by a $\pi$-circuit as shown in the Fig. 2(a) with
$\left[G_{1(2)}\right]=G_{1(2)}\left[I_{3\times 3}\right]$, where
$G_{1}=G_{\lambda}\mathrm{tanh}(t/2\lambda)$,
$G_{2}=G_{\lambda}\mathrm{csch}(t/\lambda)$, $G_{\lambda}=\sigma lw/\lambda$,
and $\left[I_{3\times 3}\right]$ is the ${3\times 3}$ identity matrix. The
spin current sources in the Fig. 2(a) in $z$-$x$-$y$ basis can be written as

$\left[I_{0,zxy}^{s}\right]=\beta
G_{0}\left(V_{1}^{c}-V_{2}^{c}\right)\left[1\;0\;0\right]^{T}$ (2)

where $G_{0}=\sigma tw/l$ and $\beta=\theta_{sh}l/t$. After converting to
$(\mathbf{m},\mathbf{m^{\prime}},\mathbf{m}\times\mathbf{m^{\prime}})$ basis,
the Equation (2) becomes (see supplementary material)

$\left[I_{0}^{s}\right]=I_{0}^{s0}m_{z}\mathbf{m}+I_{0}^{s0}\frac{1}{\omega}\mathbf{m}\times\mathbf{m^{\prime}}$
(3)

where $I_{0}^{s0}=\beta G_{0}\left(V_{1}^{c}-V_{2}^{c}\right)$.

Applying KCL at the nodes 3 and 4 in the spin-circuit shown in the Fig. 2(a),
we get (see supplementary material)

$\left[{\begin{array}[]{cc}\left[G_{1}\right]+\left[G_{2}\right]&-\left[G_{2}\right]\\\
-\left[G_{2}\right]&\left[G_{1}\right]+\left[G_{2}\right]+\left[G^{\uparrow\downarrow}\right]\\\
\end{array}}\right]\left[{\begin{array}[]{c}\left[V_{3}^{s}\right]\\\
\left[V_{4}^{s}\right]\\\
\end{array}}\right]=\left[{\begin{array}[]{c}-\left[I_{0}^{s}\right]\\\
\left[I_{0}^{s}\right]+\left[I_{sp}\right]\\\ \end{array}}\right].$ (4)

Using the superposition principle in circuit theory, i.e., the contributions
due to different voltage/current sources can be considered separately, we can
determine the contributions due to SMR and spin pumping accordingly.

Figure 2:  (a) Spin-circuit representation of ST-FMR. Spin potentials are
developed at the surfaces marked by 3 and 4. $\left[I_{0}^{s}\right]$ is the
charge current dependent spin current source due to SHE and
$\left[V_{sp}\right]$ acts as the spin battery due to spin pumping. The bulk
diffusion in the SHE layer can be represented by a $\pi$-circuit defined in
the text. (b) The charge circuit for the spin-to-charge conversion by ISHE
with the spin current dependent charge current sources $I_{0}^{c}$ with
$G_{0}$ representing the conductance of the SHE layer. (c) Spin-circuit
representation of ST-FMR without considering the spin battery
$\left[V_{sp}\right]$. (d) Spin-circuit representation of ST-FMR without
considering the spin current sources $\left[I_{0}^{s}\right]$.

Figure 2(b) shows the charge circuit counterpart comprising the spin current
dependent charge current sources due to ISHE and the conductance $G_{0}$ of
the SHE layer. A conductance in case the magnet is metallic can be included by
putting the conductance in parallel to $G_{0}$ to account for the current
shunting. Roy (2017a) The charge current sources in the Fig. 2(b) can be
expressed as

$I_{0}^{c}=\beta G_{0}\left(V_{3}^{z}-V_{4}^{z}\right).$ (5)

Applying KCL at node 1 of the charge-circuit in Fig. 2(b), we get

$I_{1}^{c}=-I_{0}^{c}+G_{0}\left(V_{1}^{c}-V_{2}^{c}\right).$ (6)

The contribution due to SMR is given by the spin-circuit in the Fig. 2(c) as
follows.

$\left[{\begin{array}[]{cc}\left[G_{1}\right]+\left[G_{2}\right]&-\left[G_{2}\right]\\\
-\left[G_{2}\right]&\left[G_{1}\right]+\left[G_{2}\right]+\left[G^{\uparrow\downarrow}\right]\\\
\end{array}}\right]\left[{\begin{array}[]{c}\left[V_{3}^{s}\right]\\\
\left[V_{4}^{s}\right]\\\
\end{array}}\right]=\left[{\begin{array}[]{c}-\left[I_{0}^{s}\right]\\\
\left[I_{0}^{s}\right]\\\ \end{array}}\right].$ (7)

After solving the above equations in the matrix form, we get

$V_{3}^{m}-V_{4}^{m}=-\frac{2\,\mathrm{tanh}(t/2\lambda)}{G_{\lambda}}I_{0}^{s0}m_{z},$
(8)
$V_{3}^{m^{\prime}}-V_{4}^{m^{\prime}}=\frac{1}{G_{\lambda}}Im\\{F(t)\\}I_{0}^{s0}\frac{1}{\omega},$
(9)

and

$\displaystyle V_{3}^{m\times m^{\prime}}-V_{4}^{m\times
m^{\prime}}=-\frac{\mathrm{tanh}(t/\lambda)}{G_{\lambda}}I_{0}^{s0}\frac{1}{\omega}$
$\displaystyle-\frac{(G_{\lambda}\mathrm{tanh}(t/\lambda)+2G_{r})\,\mathrm{tanh}^{2}(t/\lambda)\,\mathrm{tanh}^{2}(t/2\lambda)}{(G_{\lambda}\mathrm{tanh}(t/\lambda)+2G_{r})^{2}+4G_{i}^{2}}I_{0}^{s0}\frac{1}{\omega}$
(10)

where

$F(t)=\frac{(2G^{\uparrow\downarrow}/G_{\lambda})\,\mathrm{tanh}^{2}(t/2\lambda)}{1+(2G^{\uparrow\downarrow}/G_{\lambda})\,\mathrm{coth}(t/\lambda)}.$
(11)

We can write the spin voltage differences in the longitudinal and transverse
directions as follows (see supplementary material).

$V_{3}^{z}-V_{4}^{z}=(V_{3}^{m}-V_{4}^{m})m_{z}+(V_{3}^{m\times
m^{\prime}}-V_{4}^{m\times m^{\prime}})\left(1-m_{z}^{2}\right)\omega,$ (12)
$\displaystyle V_{3}^{x}-V_{4}^{x}=$
$\displaystyle(V_{3}^{m}-V_{4}^{m})m_{x}-(V_{3}^{m^{\prime}}-V_{4}^{m^{\prime}})\,m_{y}\,\omega$
$\displaystyle-(V_{3}^{m\times m^{\prime}}-V_{4}^{m\times
m^{\prime}})\,m_{z}m_{x}\,\omega.$ (13)

From Equations (6) and (12), and using the equality
$2\mathrm{tanh}(t/2\lambda)=\mathrm{tanh}(t/\lambda)\,(1+\mathrm{tanh}^{2}(t/2\lambda))$,
we get

$\displaystyle\frac{I_{1}^{c}}{tw}=\frac{\sigma}{l}\left(V_{1}^{c}-V_{2}^{c}\right)$
$\displaystyle\left\\{1+\theta_{SH}^{2}(\lambda/t)\left[2\mathrm{tanh}(t/2\lambda)\right.\right.$
$\displaystyle-\left.\left.(1-m_{z}^{2})\,Re\\{F(t)\\}\right]\right\\}.$ (14)

Using $I_{5}^{c}=\beta G_{0}\left(V_{3}^{x}-V_{4}^{x}\right)$ and Equation
(13), we get

$\frac{I_{5}^{c}}{lt}=-\theta_{SH}^{2}\frac{\sigma\lambda}{lt}\left(V_{1}^{c}-V_{2}^{c}\right)\left[m_{z}m_{x}\,Re+m_{y}\,Im\right]\,F(t).$
(15)

The Equations (14) and (15) are known as _longitudinal_ and _transverse_ SMR,
respectively and these results match the existing results in the literature.
Chen  _et al._ (2013); Hong, Sayed, and Datta (2016) Note that the SMR is
proportional to the square of the spin Hall angle ($\theta_{sh}^{2}$) due to
the coordinated action of SHE and ISHE, and therefore materials with higher
values of $\theta_{sh}$ would aid in increasing the SMR.

Figure 3:  (a) A spin-valve structure employing SMR i.e., a trilayer with an
SHE layer sandwiched by two magnetic insulator layers. The SMR will be
dependent on the magnetization directions of the two magnetic layers. (b) The
corresponding spin-circuit representation employing the two interfacial spin-
mixing conductances $G_{1}^{\uparrow\downarrow}$ and
$G_{2}^{\uparrow\downarrow}$. and two spin batteries $[V_{sp}^{1}]$ and
$[V_{sp}^{2}]$. The spin-circuit representation of the SHE layer remains same
as shown in the Fig. 2(a). Charge potentials are developed at the surfaces
marked by 1 and 2, and marked by 5 and 6 and the charge-circuit as presented
in the Fig. 2(b) remains same.

The contribution due to spin pumping is given by the spin-circuit in the Fig.
2(d) as follows.

$\left[{\begin{array}[]{cc}\left[G_{1}\right]+\left[G_{2}\right]&-\left[G_{2}\right]\\\
-\left[G_{2}\right]&\left[G_{1}\right]+\left[G_{2}\right]+\left[G^{\uparrow\downarrow}\right]\\\
\end{array}}\right]\left[{\begin{array}[]{c}\left[V_{3}^{s}\right]\\\
\left[V_{4}^{s}\right]\\\ \end{array}}\right]=\left[{\begin{array}[]{c}0\\\
\left[I_{sp}\right]\\\ \end{array}}\right].$ (16)

After solving the above equations in the matrix form, we get

$\left.V_{3}^{m^{\prime}}-V_{4}^{m^{\prime}}\right|_{sp}=(\hbar/2e)Im\\{\eta(t)\\},$
(17)

and

$\left.V_{3}^{m\times m^{\prime}}-V_{4}^{m\times
m^{\prime}}\right|_{sp}=(\hbar/2e)Re\\{\eta(t)\\}$ (18)

where

$\eta(t)=\frac{(2G^{\uparrow\downarrow}/G_{\lambda})\,\mathrm{tanh}(t/2\lambda)}{1+(2G^{\uparrow\downarrow}/G_{\lambda})\,\mathrm{coth}(t/\lambda)}.$
(19)

Note that there is no contribution in the direction of $\mathbf{m}$ since we
are not considering any spin-wave propagation or magnons. From Equation (5),
we get

$\left.\frac{I_{0}^{c}}{tw}\right|_{sp}=\frac{\hbar}{2e}\theta_{SH}\frac{\sigma}{t}\left[Re\\{\eta(t)\\}\left(\mathbf{m}\times\mathbf{m^{\prime}}\right)_{z}+Im\\{\eta(t)\\}\mathbf{m^{\prime}}_{z}\right],$
(20)

which matches the mathematical expression derived in literature. Chiba, Bauer,
and Takahashi (2014) Note that in the ST-FMR experiments, the longitudinal SMR
and spin pumping contributions can be separated from the involved angular
dependence.

For a more complex system e.g., a spin-valve employing SMR Chen  _et al._
(2013) as shown in the Fig. 3, the analytical expressions become tedious and
this clearly signifies the prowess of the spin-circuit approach. We can simply
write a netlist [see the node numbers in squares in the Fig. 3(b)] as follows
and solve that using a circuit solver.

⬇

conductances = [1 2 $G_{1}^{\uparrow\downarrow}$; 3 4
$G_{2}^{\uparrow\downarrow}$;

2 0 $[G_{1}]$; 2 3 $[G_{2}]$; 3 0 $[G_{1}]$];

voltageSources = [1 0 $[V_{sp}^{1}]$; 4 0 $[V_{sp}^{2}]$];

currentSources = [0 2 $[I_{0}^{s}]$; 3 0 $[I_{0}^{s}]$];

While using a circuit solver, there is also the other way that we can follow
i.e., convert $G_{1}^{\uparrow\downarrow}$ and $G_{2}^{\uparrow\downarrow}$ to
$z$-$x$-$y$ basis, but mathematically working with the $z$-$x$-$y$ basis is
tedious and therefore is not followed for the derivation of analytical
expressions here.

To summarize, we have developed the spin-circuit representation of ST-FMR and
have shown that both the SMR and spin pumping contributions match the
established mathematical expressions in literature. Such circuits when complex
can be simply solved using a circuit solver and therefore helps in analysis,
benchmarking experimental results, and proposing creative designs of
functional devices. We do need to consider the interfacial spin memory loss,
Roy (2020); Gupta  _et al._ (2020); Flores  _et al._ (2020); Dolui and Nikolić
(2017); Belashchenko, Kovalev, and van Schilfgaarde (2016); Rojas-Sánchez  _et
al._ (2014) which lowers the conductivity. Recent calculations Roy (2020,
2017b, 2018); Tao  _et al._ (2018); Yu  _et al._ (2018); Berger  _et al._
(2018) signify the Elliott-Yafet spin relaxation mechanism in heavy metals for
which $\lambda$ is dependent on thickness as conductivity varies with
thickness of the samples. ST-FMR can allow us to understand and estimate the
relevant parameters in the system and such understandings can benefit the
device design using SHE, Roy (2014); Ramaswamy  _et al._ (2018) which has
potential for building future spintronic devices, alongwith other promising
energy-efficient emerging devices. Roy (2016)

See the supplementary material for detailed derivation.

This work was supported by Science and Engineering Research Board (SERB) of
India via sanction order SRG/2019/002166.

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Silsbee, Janossy, and Monod (1979) R. H. Silsbee, A. Janossy, and P. Monod, Phys. Rev. B 19, 4382 (1979).
  * Heinrich  _et al._ (1987) B. Heinrich, K. B. Urquhart, A. S. Arrott, J. F. Cochran, K. Myrtle, and S. T. Purcell, Phys. Rev. Lett. 59, 1756 (1987).
  * Guan  _et al._ (2007) Y. Guan, W. E. Bailey, E. Vescovo, C.-C. Kao, and D. A. Arena, J. Magn. Magn. Mater. 312, 374 (2007).
  * Wang, Ramaswamy, and Yang (2018) Y. Wang, R. Ramaswamy, and H. Yang, J. Phys. D: Appl. Phys 51, 273002 (2018).
  * Landau and Lifshitz (1935) L. Landau and E. Lifshitz, Phys. Z. Sowjet. 8, 153 (1935).
  * Gilbert (2004) T. L. Gilbert, IEEE Trans. Magn. 40, 3443 (2004).
  * Dyakonov and Perel (1971) M. I. Dyakonov and V. I. Perel, Sov. Phys. JETP Lett. 13, 467 (1971).
  * Hirsch (1999) J. E. Hirsch, Phys. Rev. Lett. 83, 1834 (1999).
  * Zhang (2000) S. Zhang, Phys. Rev. Lett. 85, 393 (2000).
  * Murakami, Nagaosa, and Zhang (2003) S. Murakami, N. Nagaosa, and S. C. Zhang, Science 301, 1348 (2003).
  * Sinova  _et al._ (2004) J. Sinova, D. Culcer, Q. Niu, N. A. Sinitsyn, T. Jungwirth, and A. H. MacDonald, Phys. Rev. Lett. 92, 126603 (2004).
  * Sinova  _et al._ (2015) J. Sinova, S. O. Valenzuela, J. Wunderlich, C. H. Back, and T. Jungwirth, Rev. Mod. Phys. 87, 1213 (2015).
  * Valenzuela and Tinkham (2006) S. O. Valenzuela and M. Tinkham, Nature 442, 176 (2006).
  * Roy (2014) K. Roy, J. Phys. D: Appl. Phys. 47, 422001 (2014).
  * Wang  _et al._ (2006) X. Wang, G. E. W. Bauer, B. J. van Wees, A. Brataas, and Y. Tserkovnyak, Phys. Rev. Lett. 97, 216602 (2006).
  * Jiao and Bauer (2013) H. Jiao and G. E. W. Bauer, Phys. Rev. Lett. 110, 217602 (2013).
  * Azevedo  _et al._ (2005) A. Azevedo, L. H. Vilela-Leão, R. L. Rodríguez-Suárez, A. B. Oliveira, and S. M. Rezende, J. Appl. Phys. 97, 10C715 (2005).
  * Saitoh  _et al._ (2006) E. Saitoh, M. Ueda, H. Miyajima, and G. Tatara, Appl. Phys. Lett. 88, 182509 (2006).
  * Tserkovnyak, Brataas, and Bauer (2002a) Y. Tserkovnyak, A. Brataas, and G. E. W. Bauer, Phys. Rev. Lett. 88, 117601 (2002a).
  * Tserkovnyak, Brataas, and Bauer (2002b) Y. Tserkovnyak, A. Brataas, and G. E. W. Bauer, Phys. Rev. B 66, 224403 (2002b).
  * Tserkovnyak  _et al._ (2005) Y. Tserkovnyak, A. Brataas, G. E. W. Bauer, and B. I. Halperin, Rev. Mod. Phys. 77, 1375 (2005).
  * Brataas  _et al._ (2012) A. Brataas, Y. Tserkovnyak, G. E. W. Bauer, and P. J. Kelly, “Spin pumping and spin transfer,” in  _Spin Current_ , Vol. 17, edited by S. Maekawa, S. O. Valenzuela, E. Saitoh, and T. Kimura (Oxford University Press, 2012) pp. 87–135.
  * Slonczewski (1996) J. C. Slonczewski, J. Magn. Magn. Mater. 159, L1 (1996).
  * Berger (1996) L. Berger, Phys. Rev. B 54, 9353 (1996).
  * Onsager (1931a) L. Onsager, Phys. Rev. 37, 405 (1931a).
  * Onsager (1931b) L. Onsager, Phys. Rev. 38, 2265 (1931b).
  * Liu  _et al._ (2011) L. Liu, T. Moriyama, D. C. Ralph, and R. A. Buhrman, Phys. Rev. Lett. 106, 036601 (2011).
  * Wang  _et al._ (2014) Y. Wang, P. Deorani, X. Qiu, J. H. Kwon, and H. Yang, Appl. Phys. Lett. 105, 152412 (2014).
  * Zhang  _et al._ (2015) W. Zhang, W. Han, X. Jiang, S. H. Yang, and S. S. P. Parkin, Nature Phys. 11, 496 (2015).
  * Jungfleisch  _et al._ (2017) M. B. Jungfleisch, J. Ding, W.Zhang, W. Jiang, J. E. Pearson, V. Novosad, and A. Hoffmann, Nano Lett. 17, 8 (2017).
  * Chen  _et al._ (2013) Y. T. Chen, S. Takahashi, H. Nakayama, M. Althammer, S. T. B. Goennenwein, E. Saitoh, and G. E. W. Bauer, Phys. Rev. B 87, 144411 (2013).
  * Nakayama  _et al._ (2013) H. Nakayama, M. Althammer, Y. T. Chen, K. Uchida, Y. Kajiwara, D. Kikuchi, T. Ohtani, S. Geprägs, M. Opel, S. Takahashi, R. Gross, G. E. W. Bauer, S. T. B. Goennenwein, and E. Saitoh, Phys. Rev. Lett. 110, 206601 (2013).
  * Chiba, Bauer, and Takahashi (2014) T. Chiba, G. E. W. Bauer, and S. Takahashi, Phys. Rev. Appl. 2, 034003 (2014).
  * Chiba  _et al._ (2015) T. Chiba, M. Schreier, G. E. W. Bauer, and S. Takahashi, J. Appl. Phys. 117, 17C715 (2015).
  * Schreier  _et al._ (2015) M. Schreier, T. Chiba, A. Niedermayr, J. Lotze, H. Huebl, S. Geprägs, S. Takahashi, G. E. W. Bauer, R. Gross, and S. T. B. Goennenwein, Phys. Rev. B 92, 144411 (2015).
  * Jungfleisch  _et al._ (2015) M. B. Jungfleisch, A. V. Chumak, A. Kehlberger, V. Lauer, D. H. Kim, M. C. Onbasli, C. A. Ross, M. Kläui, and B. Hillebrands, Phys. Rev. B 91, 134407 (2015).
  * Sklenar  _et al._ (2015) J. Sklenar, W. Zhang, M. B. Jungfleisch, W. Jiang, H. Chang, J. E. Pearson, M. Wu, J. B. Ketterson, and A. Hoffmann, Phys. Rev. B 92, 174406 (2015).
  * Thomson (1856) W. Thomson, Philos. Trans. R. Soc. London 146, 649 (1856).
  * Baibich  _et al._ (1988) M. N. Baibich, J. M. Broto, A. Fert, F. NguyenVanDau, F. Petroff, P. Etienne, G. Creuzet, A. Friederich, and J. Chazelas, Phys. Rev. Lett. 61, 2472 (1988).
  * Binasch  _et al._ (1989) G. Binasch, P. Grünberg, F. Saurenbach, and W. Zinn, Phys. Rev. B 39, 4828 (1989).
  * Julliere (1975) M. Julliere, Phys. Lett. A 54, 225 (1975).
  * Butler  _et al._ (2001) W. H. Butler, X. G. Zhang, T. C. Schulthess, and J. M. MacLaren, Phys. Rev. B 63, 054416 (2001).
  * Parkin  _et al._ (2004) S. S. P. Parkin, C. Kaiser, A. Panchula, P. M. Rice, B. Hughes, M. G. Samant, and S. H. Yang, Nature Mater. 3, 862 (2004).
  * Yuasa  _et al._ (2004) S. Yuasa, T. Nagahama, A. Fukushima, Y. Suzuki, and K. Ando, Nature Mater. 3, 868 (2004).
  * Rabaey, Chandrakasan, and Nikoliç (2003) J. M. Rabaey, A. P. Chandrakasan, and B. Nikoliç,  _Digital Integrated Circuits_ (Pearson Education, 2003).
  * (46) _HSPICE, Circuit simulation software_ , Synopsys, www.synopsys.com.
  * Sah (1971) C.-T. Sah, Phys. Stat. Sol. (a) 7, 541 (1971).
  * Brataas, Bauer, and Kelly (2006) A. Brataas, G. E. W. Bauer, and P. J. Kelly, Phys. Rep. 427, 157 (2006).
  * Srinivasan  _et al._ (2014) S. Srinivasan, V. Diep, B. Behin-Aein, A. Sarkar, and S. Datta, “Modeling multi-magnet networks interacting via spin currents,” in  _Handbook of Spintronics_, edited by Y. Xu, D. D. Awschalom, and J. Nitta (Springer Netherlands, Dordrecht, 2014) pp. 1–49.
  * Hong, Sayed, and Datta (2016) S. Hong, S. Sayed, and S. Datta, IEEE Trans. Nanotech. 15, 225 (2016).
  * Roy (2017a) K. Roy, Phys. Rev. Appl. (Lett.) 8, 011001 (2017a).
  * Brouwer (1998) P. W. Brouwer, Phys. Rev. B 58, R10135 (1998).
  * Büttiker, Thomas, and Prêtre (1994) M. Büttiker, H. Thomas, and A. Prêtre, Z. Phys. B 94, 133 (1994).
  * Roy (2020) K. Roy, Appl. Phys. Lett. 117, 022404 (2020).
  * Gupta  _et al._ (2020) K. Gupta, R. J. H. Wesselink, R. Liu, Z. Yuan, and P. J. Kelly, Phys. Rev. Lett. 124, 087702 (2020).
  * Flores  _et al._ (2020) G. G. B. Flores, A. A. Kovalev, M. van Schilfgaarde, and K. D. Belashchenko, Phys. Rev. B 101, 224405 (2020).
  * Dolui and Nikolić (2017) K. Dolui and B. K. Nikolić, Phys. Rev. B 96, 220403 (2017).
  * Belashchenko, Kovalev, and van Schilfgaarde (2016) K. D. Belashchenko, A. A. Kovalev, and M. van Schilfgaarde, Phys. Rev. Lett. 117, 207204 (2016).
  * Rojas-Sánchez  _et al._ (2014) J. C. Rojas-Sánchez, N. Reyren, P. Laczkowski, W. Savero, J. P. Attané, C. Deranlot, M. Jamet, J. M. George, L. Vila, and H. Jaffrès, Phys. Rev. Lett. 112, 106602 (2014).
  * Roy (2017b) K. Roy, Phys. Rev. B 96, 174432 (2017b).
  * Roy (2018) K. Roy, Proc. SPIE Nanoscience (Spintronics XI) 10732, 1073207 (2018).
  * Tao  _et al._ (2018) X. Tao, Q. Liu, B. Miao, R. Yu, Z. Feng, L. Sun, B. You, J. Du, K. Chen, S. Zhang, L. Zhang, Z. Yuan, D. Wu, and H. Ding, Sci. Adv. 4, eaat1670 (2018).
  * Yu  _et al._ (2018) R. Yu, B. F. Miao, L. Sun, Q. Liu, J. Du, P. Omelchenko, B. Heinrich, M. Wu, and H. F. Ding, Phys. Rev. Mater. 2, 074406 (2018).
  * Berger  _et al._ (2018) A. J. Berger, E. R. J. Edwards, H. T. Nembach, O. Karis, M. Weiler, and T. J. Silva, Phys. Rev. B 98, 024402 (2018).
  * Ramaswamy  _et al._ (2018) R. Ramaswamy, J. M. Lee, K. Cai, and H. Yang, Appl. Phys. Rev. 5, 031107 (2018).
  * Roy (2016) K. Roy, SPIN 6, 1630001 (2016).


# Structural relaxation and delayed yielding in cyclically sheared Cu-Zr
metallic glasses

Nikolai V. Priezjev1,2 1Department of Civil and Environmental Engineering,
Howard University, Washington, D.C. 20059  2Department of Mechanical and
Materials Engineering, Wright State University, Dayton, OH 45435

###### Abstract

The yielding transition, structural relaxation, and mechanical properties of
metallic glasses subjected to repeated loading are examined using molecular
dynamics simulations. We consider a poorly-annealed Cu-Zr amorphous alloy
periodically deformed in a wide range of strain amplitudes at room
temperature. It is found that low-amplitude cyclic loading leads to a
logarithmic decay of the potential energy, and lower energy states are
attained when the strain amplitude approaches a critical point from below.
Moreover, the potential energy after several thousand loading cycles is a
linear function of the peak value of the stress overshoot during startup
continuous shear deformation of the annealed sample. We show that the process
of structural relaxation involves collective, irreversible rearrangements of
groups of atoms whose spatial extent is most pronounced at the initial stage
of loading and higher strain amplitudes. At the critical amplitude, the glass
becomes mechanically annealed for a number of transient cycles and then yields
via formation of a shear band. The yielding transition is clearly marked by
abrupt changes in the potential energy, storage modulus, and fraction of atoms
with large nonaffine displacements.

Keywords: yielding transition, metallic glasses, plastic deformation, cyclic
loading, molecular dynamics simulations

##  I Introduction

Establishing structure-property-performance relations for bulk metallic
glasses is important for various structural and functional applications Gao22
; Jiang21 . Owing to their amorphous atomic structure, metallic glasses offer
a number of unique properties, such as high strength, large elastic limit, as
well as superior wear and corrosion resistance Egami13 ; Qiao19 . An
outstanding challenge that limits their widespread use is that sufficiently
well annealed glasses are prone to brittle failure via formation of nanoscale
shear bands Hufnagel16 . In addition, metallic glasses subjected cyclic
variations in stresses or strains exhibit relatively low fatigue limit and
fatigue life, which can be affected by the sample size, chemical composition,
cycling frequency, and surface conditions Liaw18 ; Menzel06 ; Sha2020 . To
enhance their ductility, metallic glasses can be rejuvenated using several
thermo-mechanical processing methods, including high-pressure torsion, ion
irradiation, thermal cycling, cold rolling, and elastostatic loading Greer16 ;
Ketov15 . Alternatively, an atomic structure can be simply reset by heating
above the glass transition temperature and subsequent rapid cooling to the
glass state or a supercooled glass former can be frozen under applied stress
Mota21 ; PriezMETAL21 . In recent years, laser-based additive manufacturing
techniques were also introduced to fabricate large-scale complex structures
and patient-specific implants for biomedical applications Reed23 . Despite
these advances, however, development of novel processing methods to rejuvenate
metallic glasses and improve their mechanical properties remains a challenging
task.

In the last decade, a number of molecular dynamics (MD) simulation studies
have investigated the yielding behavior and structural relaxation in amorphous
solids subjected to oscillatory deformation Priezjev13 ; Reichhardt13 ;
Sastry13 ; IdoNature15 ; Shi15 ; GaoNano15 ; Priezjev16 ; Kawasaki16 ;
Priezjev16a ; Sastry17 ; Priezjev17 ; Hecke17 ; Priezjev18 ; Priezjev18a ;
Sastry19band ; NVP18strload ; PriezSHALT19 ; Priez20ba ; KawBer20 ; NVP20altY
; Priez20delay ; BhaSastry21 ; Priez21var ; PriezCMS21 ; Peng22 ; Sastry22 ;
PriezJNCS22 ; Barrat23 ; PriezCMS23 ; Karmakar23 ; PriezMet23 ; Egami24 ;
Tang24 . Notably, it was demonstrated that in the athermal limit, binary
Lennard-Jones (LJ) glasses under low-amplitude loading gradually evolve into
periodic limit cycles with exactly reversible trajectories of atoms
Reichhardt13 ; IdoNature15 , and the potential energy of such states is lower
for better annealed glasses and higher strain amplitudes Sastry13 ; Sastry17 .
By contrast, periodic shear deformation in combination with thermal noise
facilitate collective, irreversible rearrangements of atoms and prolonged
structural relaxation Sastry22 ; PriezJNCS22 . In addition, lower energy
states can be accessed when cyclic shear is periodically alternated along two
or three mutually perpendicular planes PriezSHALT19 ; Karmakar23 . Moreover,
the range of energy states attainable in thermal glasses during low-amplitude
loading can be extended by increasing strain amplitude slightly above a
critical point every few cycles Priez21var . Furthermore, when the loading
amplitude exceeds a critical value, amorphous alloys undergo a yielding
transition after a number of transient cycles GaoNano15 ; Sastry17 ;
Priezjev17 ; Priezjev18a ; Sastry19band ; Priez20ba ; NVP20altY ; Priez20delay
; PriezCMS21 ; PriezJNCS22 ; PriezCMS23 ; PriezMet23 . Interestingly, the
critical strain amplitude in athermal systems might depend on the glass
stability KawBer20 ; BhaSastry21 , whereas at about half $T_{g}$, the yielding
transition is delayed in better annealed glasses but the critical strain
amplitude remains unchanged PriezJNCS22 . It was recently shown that well
annealed binary glasses fail via shear band formation, and number of cycles to
reach the yielding transition is well described by the power-law function of
the difference between the strain amplitude and its critical value PriezCMS23
; PriezMet23 . However, the role of loading conditions and preparation history
on the critical behavior of metallic glasses is not yet fully understood.

In this paper, the yielding behavior and structural relaxation in a Cu-Zr
metallic glass under periodic shear deformation are studied using molecular
dynamics simulations. We consider the binary glass that is first rapidly
cooled across the glass transition to room temperature and then cyclically
loaded in a wide range of strain amplitudes around a critical value. It will
be shown that low-amplitude loading leads to a logarithmic decay of the
potential energy during thousands of cycles. When the loading amplitude
increases toward the critical value, the average size of plastically deformed
domains becomes larger and the glass is relocated to lower energy states. On
the contrary, we find that the yielding transition at the critical strain
amplitude is marked by abrupt changes in the storage modulus, potential
energy, and the number of atoms with large nonaffine displacements, which are
localized within a narrow shear band.

The rest of this paper is structured as follows. The MD simulation setup,
parameter values, as well as cooling and deformation protocols are described
in the next section. The analysis of the potential energy, shear stress, and
nonaffine displacements are presented in section III. A summary of the results
is provided in the last section.

##  II MD simulations

In our study, the metallic glass was represented by a binary mixture of Cu and
Zr atoms, which interacted via the embedded atom method (EAM) potentials
CuZrEAM09 ; Ma11 . The total number of atoms in the
$\text{Cu}_{50}\text{Zr}_{50}$ glass is $60\,000$. The preparation procedure
consisted of several steps. First, the binary mixture was thoroughly
equilibrated in a periodic box at the temperature of $2000\,\text{K}$ and zero
pressure. This temperature is well above the glass transition temperature
$T_{g}\approx 675\,\text{K}$. We followed the cooling protocol by Fan and Ma
FanMa21 , where the Cu-Zr system is initially cooled to $1500\,\text{K}$ at
the rate of $10^{13}\,\text{K/s}$ and then to $300\,\text{K}$ at
$10^{12}\,\text{K/s}$. As a result, the effective cooling rate was
$10^{12}\,\text{K/s}$. The sample was prepared using the Nosé-Hoover
thermostat, zero applied pressure, and periodic boundary conditions. The
equations of motion were integrated with the time step $\triangle
t=1.0\,\text{fs}$ Lammps .

The periodic shear deformation was applied along the $xz$ plane as follows:

$\gamma_{xz}(t)=\gamma_{0}\,\text{sin}(2\pi t/T),$ (1)

where $\gamma_{0}$ is the strain amplitude and $T$ is the oscillation period.
In our setup, the oscillation period was set to $1.0\,\text{ns}$ ($10^{6}$ MD
time steps), and the strain amplitude was varied in the range
$0.020\leqslant\gamma_{0}\leqslant 0.061$. The MD simulations were carried out
in the NVT ensemble with the temperature of $300\,\text{K}$ and the linear
size of the periodic box of $101.8\,\text{\AA}$. A typical simulation of 4000
shear cycles took about 6000 hours using 40 processors in parallel. Due to
computational limitations, the data for the potential energy, shear stress,
and atomic configurations were collected for only one independent sample.

##  III Results

One of the key factors that strongly affects mechanical properties of metallic
glasses is the rate of cooling during glass formation Greer16 . In general,
glasses obtained by cooling at a slower rate become more stable and, upon
deformation, exhibit a stress overshoot followed by plastic flow, whereas
rapidly cooled glasses are settled at higher energy states and yield more
smoothly Egami13 . Alternatively, it was recently demonstrated that poorly
annealed LJ glasses can be relocated to lower energy states via low-amplitude
cyclic loading at a finite temperature below $T_{g}$ Sastry22 ; PriezJNCS22 .
However, the resulting change in the potential energy at given strain
amplitude and temperature as well as the yielding behavior near the critical
strain amplitude for more realistic models of glasses remain to be determined.
In the present study, we considered the $\text{Cu}_{50}\text{Zr}_{50}$ glass
that was first rapidly cooled to room temperature and then periodically
strained for $4000$ shear cycles over a broad range of strain amplitudes. This
relatively large number of loading cycles was set based on the results of the
previous MD study where a rapidly cooled binary LJ glass (60 000 atoms)
yielded only after about $2500$ cycles at a critical strain amplitude
NVP20altY .

We first plot potential energy minima after each cycle at zero strain in Fig.
1 for strain amplitudes $0.020\leqslant\gamma_{0}\leqslant 0.055$. As is
evident, the glass is relocated to progressively lower energy states upon
continued loading, and cycle-to-cycle fluctuations are enhanced at larger
strain amplitudes. Moreover, cyclic loading at higher strain amplitudes (up to
a critical value) allows for rearrangement of larger clusters of atoms during
each cycle, leading, on average, to more relaxed states. These conclusions are
in agreement with the results of previous MD studies of binary LJ glasses
Sastry17 ; Priezjev18 ; Sastry22 ; PriezJNCS22 . In Fig. 2, the same data for
the potential energy as a function of the cycle number are replotted on a
semi-log scale. It can be clearly observed that for each strain amplitude, the
potential energy closely follows a logarithmic decay when $t/T\gtrsim 10$,
suggesting a possibility of reaching lower energy states upon further loading.
Interestingly, periodic deformation for 4000 shear cycles at
$\gamma_{0}=0.055$ resulted in the potential energy
$U\approx-4.951\,\text{eV}$ (see Fig. 2), which is nearly the same as in the
case of a better annealed $\text{Cu}_{50}\text{Zr}_{50}$ glass prepared at a
100 times slower cooling rate PriezMet23 . In other words, the rapidly cooled
glass ($10^{12}\,\text{K/s}$) was mechanically annealed to the potential
energy level of a more slowly cooled glass ($10^{10}\,\text{K/s}$), which was
considered in the previous MD study PriezMet23 .

The variation of the potential energy during cyclic loading at higher strain
amplitudes, $\gamma_{0}\geqslant 0.055$, is presented in Fig. 3. For
reference, the same data for $\gamma_{0}=0.055$ as in Fig. 1 are also included
in Fig. 3 for the first 1300 cycles. It can be readily seen that, following a
number of transient cycles, the potential energy abruptly increases,
indicating shear band formation at the yielding transition. Note that except
for a relatively large strain amplitude $\gamma_{0}=0.061$, rapidly cooled
glass is mechanically annealed for a number of cycles before yielding. The
number of cycles to reach the yielding transition generally increases upon
reducing strain amplitude toward a critical value. As shown in Fig. 3, the
maximum number of transient cycles is about 650 at the critical strain
amplitude $\gamma_{0}=0.056$. The transient behavior can be rationalized as
follows. When a rapidly cooled glass is initially strained up to
$\gamma_{xz}=\pm\gamma_{0}$ in Eq. (1), the maximum stress remains relatively
low and plastic deformation is homogeneously distributed within the sample.
Upon further loading, the glass becomes more stable (lower potential energy at
zero strain), the maximum stress increases, and the formation of a large-scale
plastic event becomes more probable. We finally comment that a better annealed
$\text{Cu}_{50}\text{Zr}_{50}$ glass under cyclic loading did not yield for
700 cycles at $\gamma_{0}=0.056$, when the potential energy was
$U\approx-4.951\,\text{eV}$ PriezMet23 . Instead, the critical strain
amplitude was found to be $\gamma_{0}=0.057$ for similar loading conditions
PriezMet23 .

Along with the potential energy, the time dependence of shear stress,
$\sigma_{xz}(t)$, was analyzed for different strain amplitudes. For each shear
cycle, we computed the storage modulus,
$G^{\prime}=\sigma^{max}_{xz}/\gamma_{0}\,\text{cos}(\delta)$, where $\delta$
is the phase lag between stress and strain McKinley11 . The variation of the
storage modulus is presented in Fig. 4 for strain amplitudes below the
critical value, i.e., $\gamma_{0}\leqslant 0.055$. It can be seen that for
each strain amplitude, the storage modulus increases roughly logarithmically
as a function of the cycle number, which is consistent with the decay of the
potential energy reported in Fig. 2. The largest increase in $G^{\prime}$
during 4000 cycles is found for the strain amplitude $\gamma_{0}=0.055$, which
is just below the critical value (see Fig. 4). Similar to results for binary
LJ glasses Priezjev18a , the storage modulus is larger for cyclic loading at
smaller strain amplitudes, where deviation from the elastic regime of
deformation is reduced.

We next examine the stress-strain response of the binary glass subjected to
startup continuous shear deformation with a constant strain rate of
$10^{-5}\,\text{ps}^{-1}$. The dependence of shear stress as a function of
strain is shown in Fig. 5 after loading for 4000 cycles at the indicated
strain amplitudes. For comparison, the data for the glass after rapid cooling
but before cyclic deformation are also included in Fig. 5. As expected, the
rapidly cooled glass under steady strain exhibits a smooth crossover to
plastic flow (the violet curve in Fig. 5). By contrast, plastic flow of
mechanically annealed glasses is preceded by the yielding peak, which becomes
more pronounced in glasses previously loaded at larger strain amplitudes. In
Fig. 6, we summarize results for the shear modulus, computed from the slope of
stress-strain curves at $\gamma_{xz}\leqslant 0.01$, and the peak value of the
stress overshoot as functions of the strain amplitude and the potential energy
after 4000 cycles. It can be concluded from Fig. 6 (b) that the potential
energy of the mechanically annealed glass is approximately linearly related to
the magnitude of the yielding peak. On the other hand, the shear modulus
remains rather insensitive to the loading amplitude, as shown Fig. 6 (c).

The storage modulus as a function of the cycle number is plotted in Fig. 7 for
$\gamma_{0}\geqslant 0.055$. Note that the data for the strain amplitude
$\gamma_{0}=0.055$ are the same as in Fig. 4. It is seen in Fig. 7 that in the
range $0.056\leqslant\gamma_{0}\leqslant 0.060$, cyclic loading leads to an
increase in $G^{\prime}$ for a number of transient cycles, followed by the
yielding transition, which is marked by an abrupt drop in storage modulus.
Interestingly, the storage modulus as a function of the cycle number at the
critical strain amplitude, $\gamma_{0}=0.056$, closely follows the data for
$\gamma_{0}=0.055$ during the first 600 cycles. Thus, the onset of yielding
cannot be determined from the stress variation during the initial stage of
deformation near the critical strain amplitude. As shown in Fig. 7, the number
of transient cycles is reduced at higher strain amplitudes,
$\gamma_{0}>0.056$. Overall, these results correlate well with dependence of
the potential energy minima on the number of cycles reported in Fig. 3.

At the microscopic level, plastic deformation in disordered materials can be
described via the so-called nonaffine displacements of atoms with respect to
their neighbors Falk98 . In practice, the nonaffine quantity for an atom
displaced from the position vector $\mathbf{r}_{i}(t)$ to
$\mathbf{r}_{i}(t+\Delta t)$ during the time interval $\Delta t$ is computed
by minimizing the following expression:

$D^{2}(t,\Delta
t)=\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}\Big{\\{}\mathbf{r}_{j}(t+\Delta
t)-\mathbf{r}_{i}(t+\Delta
t)-\mathbf{J}_{i}\big{[}\mathbf{r}_{j}(t)-\mathbf{r}_{i}(t)\big{]}\Big{\\}}^{2},$
(2)

where $\mathbf{J}_{i}$ is the transformation matrix and the sum is taken over
neighboring atoms that are initially located within $4.0\,\text{\AA}$ from
$\mathbf{r}_{i}(t)$. It was recently shown that in glasses under periodic
deformation, rearrangements of atoms become irreversible when their nonaffine
displacements are greater than a cage size Priezjev16 . In addition, previous
MD simulation study has indicated that the cage size is about
$0.6\,\text{\AA}$ for $\text{Cu}_{50}\text{Zr}_{50}$ metallic alloy near the
glass transition temperature Douglas19 .

The dependence of the nonaffine measure as a function of the cycle number is
presented in Fig. 8 for $\gamma_{0}\leqslant 0.055$ and in Fig. 9 for
$\gamma_{0}\geqslant 0.055$. In our analysis, the nonaffine measure was
computed at the beginning and end of each cycle, $\Delta t=T$ in Eq. (2), and
then averaged over all atoms. For cyclic loading below the critical strain
amplitude, $\gamma_{0}=0.056$, the glass evolves toward lower energy states
via a sequence of plastic events, whose size is reduced at lower amplitudes,
and, as a result, the average of $D^{2}$ is smaller at lower $\gamma_{0}$, as
shown in Fig. 8. Note that the largest fluctuations in $D^{2}$ occur at the
strain amplitude $\gamma_{0}=0.055$, whereas the lower bound, $D^{2}\approx
0.1\,\text{\AA}^{2}$, is dominated by displacements of atoms within their
cages. Further, the results in the inset to Fig. 8 for the fraction of atoms
with large nonaffine displacements demonstrate clearly that the number of
atoms involved in plastic events decays with increasing number of cycles. By
contrast, the average of $D^{2}$ increases sharply at the yielding transition
for $\gamma_{0}\geqslant 0.056$ (see Fig. 9), which is in agreement with the
critical behavior of $U$ and $G^{\prime}$ reported in Figures 3 and 7,
respectively. Moreover, after the yielding transition, plastic flow is
localized within a shear band that contains a large fraction of atoms,
$n_{f}\gtrsim 0.25$ (see Fig. 10).

The spatial organization of plastically deformed domains during the relaxation
process at the strain amplitude $\gamma_{0}=0.055$ is illustrated in Fig. 11.
The snapshots include configurations of atoms after $n$ cycles at zero strain,
and the colorcode denotes the magnitude of the nonaffine measure for $\Delta
t=T$ in Eq. (2). Note that atoms with relatively small nonaffine
displacements, $D^{2}[(n-1)\,T,T]<0.49\,\text{\AA}^{2}$, were omitted for
clarity, and, therefore, the void space represents elastically deformed
regions. As shown in Fig. 11 (a), the transition to lower energy states
initially proceeds via large irreversible displacements of atoms that form
several compact clusters homogeneously distributed in the sample. Upon further
loading, the total number of atoms with large nonaffine displacements is
significantly reduced, see Fig. 11 (b), and eventually plastic rearrangements
remain localized only in several isolated clusters, see Fig. 11 (c, d). Hence,
these results demonstrate that poorly annealed glasses subjected to prolonged
low-amplitude periodic deformation become more stable and reversible. On the
contrary, cyclic loading at the critical strain amplitude, $\gamma_{0}=0.056$,
leads to a delayed yielding transition and formation of a shear band across
the whole sample, as shown in Fig. 12. A close comparison of the results for
$\gamma_{0}=0.056$ in Figs. 10 and 12 reveals details of the yielding
transition; namely, shear band initiation [$n_{f}=0.07$ in Fig. 12 (b)],
propagation of plastic regions along the $yz$ plane [$n_{f}=0.20$ in Fig. 12
(c)], and formation of a fully developed shear band [$n_{f}=0.28$ in Fig. 12
(d)]. Together, these findings emphasize the role of plastic rearrangements in
structural relaxation and yielding of metallic glasses under cyclic loading.

##  IV Conclusions

In summary, we investigated the critical behavior and mechanical annealing of
a Cu-Zr metallic glass subjected to periodic shear deformation by means of
molecular dynamics simulations. The glass was prepared via rapid cooling from
the liquid state to room temperature and then subjected to oscillatory shear
at strain amplitudes ranged from slightly above to well below a critical
value. It was found that low-amplitude cyclic loading gradually relocates the
glass to lower energy states, and the storage modulus approximately follows a
logarithmic dependence on the number of shear cycles. In addition, we showed
that the potential energy after four thousand shear cycles depends linearly on
the yielding stress of the annealed glass under monotonic shear deformation.
The structural relaxation proceeds via a sequence of plastic rearrangements of
clusters of atoms whose typical size is larger at higher strain amplitudes. By
contrast, when the strain amplitude is greater than a critical value, the
glass is mechanically annealed for a number of cycles followed by an abrupt
increase in the potential energy due to flow localization within a shear band.
The formation of the shear band is accompanied with the increase of the
fraction of atoms with large nonaffine displacements and a drop in storage
modulus.

## Acknowledgments

The financial support from the National Science Foundation (CNS-1531923) is
gratefully acknowledged. MD simulations were carried out at the Wright State
University’s Computing Facility and the Ohio Supercomputer Center using the
LAMMPS code Lammps .

Figure 1: The potential energy at the end of each cycle versus cycle number
for the indicated strain amplitudes $\gamma_{0}$. The period of oscillation is
$T=1.0\,\text{ns}$. Figure 2: A linear-log plot of the same data as in Fig. 1.
Figure 3: The variation of potential energy when strain is zero as a function
of the number of cycles for strain amplitudes $\gamma_{0}=0.055$, $0.056$,
$0.057$, $0.060$ and $0.061$. The data for $\gamma_{0}=0.055$ are the same as
in Fig. 1.  Figure 4: The storage modulus $G^{\prime}$ (in units of GPa) as a
function of the cycle number for the indicated values of the strain amplitude
$\gamma_{0}$. The period of oscillation is $T=1.0\,\text{ns}$. Figure 5: Shear
stress, $\sigma_{xz}$ (in units of GPa), versus strain for glasses deformed at
a constant strain rate of $10^{-5}\,\text{ps}^{-1}$ after 4000 cycles at
indicated $\gamma_{0}$. The lower violet curve is the data before cyclic
loading. Figure 6: The shear modulus, $G$ (in GPa), and the peak value of the
stress overshoot, $\sigma_{Y}$ (in GPa), as functions of the strain amplitude,
$\gamma_{0}$, and the potential energy after 4000 loading cycles, $U$. The
horizontal dashed lines indicate $G$ and $\sigma_{Y}$ for the steadily sheared
glass before cyclic loading was applied. The straight red line in the panel
(b) is the best fit to the data. Figure 7: The storage modulus $G^{\prime}$
(in units of GPa) versus cycle number for strain amplitudes
$\gamma_{0}\geqslant 0.055$. The data for $G^{\prime}$ at $\gamma_{0}=0.055$
are the same as in Fig. 4. The oscillation period is $T=1.0\,\text{ns}$.
Figure 8: The average of the nonaffine quantity, $D^{2}[(n-1)\,T,T]$ (in units
of $\text{\AA}^{2}$), as a function of the number of cycles for the indicated
strain amplitudes. The oscillation period is $T=1.0\,\text{ns}$. The inset
shows the fraction of atoms with $D^{2}[(n-1)\,T,T]>0.49\,\text{\AA}^{2}$
versus cycle number for the same strain amplitudes.  Figure 9: The average of
$D^{2}[(n-1)\,T,T]$ (in units of $\text{\AA}^{2}$) versus number of cycles for
strain amplitudes $\gamma_{0}\geqslant 0.055$. The data for the strain
amplitude $\gamma_{0}=0.055$ (the black curve) are the same as in Fig. 8.
Figure 10: The fraction of atoms with large nonaffine displacements,
$D^{2}[(n-1)\,T,T]>0.49\,\text{\AA}^{2}$, as a function of the number of
cycles for the tabulated strain amplitudes. The oscillation period is
$T=1.0\,\text{ns}$. Figure 11: Selected configurations of atoms in the
$\text{Cu}_{50}\text{Zr}_{50}$ glass loaded at the strain amplitude
$\gamma_{0}=0.055$. The nonaffine quantity is (a)
$D^{2}(50\,T,T)>0.49\,\text{\AA}^{2}$, (b)
$D^{2}(500\,T,T)>0.49\,\text{\AA}^{2}$, (c)
$D^{2}(1000\,T,T)>0.49\,\text{\AA}^{2}$, and (d)
$D^{2}(3000\,T,T)>0.49\,\text{\AA}^{2}$. Cu and Zr atoms are not drawn to
scale. The legend in the panel (d) indicates the magnitude of the nonaffine
quantity.  Figure 12: Snapshots of atomic configurations of the binary glass
periodically deformed at the critical strain amplitude $\gamma_{0}=0.056$. The
nonaffine measure is (a) $D^{2}(300\,T,T)>0.49\,\text{\AA}^{2}$, (b)
$D^{2}(600\,T,T)>0.49\,\text{\AA}^{2}$, (c)
$D^{2}(660\,T,T)>0.49\,\text{\AA}^{2}$, and (d)
$D^{2}(800\,T,T)>0.49\,\text{\AA}^{2}$. Atoms are not depicted to scale. The
colorcode is the same as in Fig. 11.

## References

  * (1) K. Gao, X. G. Zhu, L. Chen, W. H. Li, X. Xu, B. T. Pan, W. R. Li, W. H. Zhou, L. Li, W. Huang, and Y. Li, Recent development in the application of bulk metallic glasses, J. Mater. Sci. Technol. 131, 115 (2022).
  * (2) H. Jiang, T. Shang, H. Xian, B. Sun, Q. Zhang, Q. Yu, H. Bai, L. Gu, and W. Wang, Functional properties of amorphous alloys, Small Struct. 2, 2000057 (2021).
  * (3) T. Egami, T. Iwashita, and W. Dmowski, Mechanical properties of metallic glasses, Metals 3, 77 (2013).
  * (4) J. C. Qiao, Q. Wang, J. M. Pelletier, H. Kato, R. Casalini, D. Crespo, E. Pineda, Y. Yao, Y. Yang, Structural heterogeneities and mechanical behavior of amorphous alloys, Prog. Mater. Sci. 104, 250 (2019).
  * (5) T. C. Hufnagel, C. A. Schuh, and M. L. Falk, Deformation of metallic glasses: Recent developments in theory, simulations, and experiments, Acta Materialia 109, 375 (2016).
  * (6) H. Jia, G. Wang, S. Chen, Y. Gao, W. Li, and P. K. Liaw, Fatigue and fracture behavior of bulk metallic glasses and their composites, Prog. Mater. Sci. 98, 168 (2018).
  * (7) B. C. Menzel and R. H. Dauskardt, Stress-life fatigue behavior of a Zr-based bulk metallic glass, Acta Mater. 54, 935 (2006).
  * (8) Z. Sha, W. Lin, L. H. Poh, G. Xing, Z. Liu, T. Wang, and H. Gao, Fatigue of metallic glasses, Appl. Mech. Rev. 72, 050801 (2020).
  * (9) Y. Sun, A. Concustell, and A. L. Greer, Thermomechanical processing of metallic glasses: Extending the range of the glassy state, Nat. Rev. Mater. 1, 16039 (2016).
  * (10) S. V. Ketov, Y. H. Sun, S. Nachum, Z. Lu, A. Checchi, A. R. Beraldin, H. Y. Bai, W. H. Wang, D. V. Louzguine-Luzgin, M. A. Carpenter, and A. L. Greer, Rejuvenation of metallic glasses by non-affine thermal strain, Nature 524, 200 (2015).
  * (11) R. M. O. Mota, E. T. Lund, S. Sohn, D. J. Browne, D. C. Hofmann, S. Curtarolo, A. van de Walle, and J. Schroers, Enhancing ductility in bulk metallic glasses by straining during cooling, Commun. Mater. 2, 23 (2021).
  * (12) N. V. Priezjev, Cooling under applied stress rejuvenates amorphous alloys and enhances their ductility, Metals 11, 67 (2021).
  * (13) A. A. Abdu Aliyu, C. Panwisawas, J. Shinjo, C. Puncreobutr, R. C. Reed, K. Poungsiri, and B. Lohwongwatana, Laser-based additive manufacturing of bulk metallic glasses: recent advances and future perspectives for biomedical applications, J. Mater. Res. Technol. 23, 2956 (2023).
  * (14) N. V. Priezjev, Heterogeneous relaxation dynamics in amorphous materials under cyclic loading, Phys. Rev. E 87, 052302 (2013).
  * (15) I. Regev, T. Lookman, and C. Reichhardt, Onset of irreversibility and chaos in amorphous solids under periodic shear, Phys. Rev. E 88, 062401 (2013).
  * (16) D. Fiocco, G. Foffi, and S. Sastry, Oscillatory athermal quasistatic deformation of a model glass, Phys. Rev. E 88, 020301(R) (2013).
  * (17) I. Regev, J. Weber, C. Reichhardt, K. A. Dahmen, and T. Lookman, Reversibility and criticality in amorphous solids, Nat. Commun. 6, 8805 (2015).
  * (18) J. Luo, K. Dahmen, P. K. Liaw, and Y. Shi, Low-cycle fatigue of metallic glass nanowires, Acta Mater. 87, 225 (2015).
  * (19) Z. D. Sha, S. X. Qu, Z. S. Liu, T. J. Wang, and H. Gao, Cyclic deformation in metallic glasses, Nano Lett. 15, 7010 (2015).
  * (20) N. V. Priezjev, Reversible plastic events during oscillatory deformation of amorphous solids, Phys. Rev. E 93, 013001 (2016).
  * (21) T. Kawasaki and L. Berthier, Macroscopic yielding in jammed solids is accompanied by a non-equilibrium first-order transition in particle trajectories, Phys. Rev. E 94, 022615 (2016).
  * (22) N. V. Priezjev, Nonaffine rearrangements of atoms in deformed and quiescent binary glasses, Phys. Rev. E 94, 023004 (2016).
  * (23) P. Leishangthem, A. D. S. Parmar, and S. Sastry, The yielding transition in amorphous solids under oscillatory shear deformation, Nat. Commun. 8, 14653 (2017).
  * (24) N. V. Priezjev, Collective nonaffine displacements in amorphous materials during large-amplitude oscillatory shear, Phys. Rev. E 95, 023002 (2017).
  * (25) S. Dagois-Bohy, E. Somfai, B. P. Tighe, and M. van Hecke, Softening and yielding of soft glassy materials, Soft Matter 13, 9036 (2017).
  * (26) N. V. Priezjev, Molecular dynamics simulations of the mechanical annealing process in metallic glasses: Effects of strain amplitude and temperature, J. Non-Cryst. Solids 479, 42 (2018).
  * (27) N. V. Priezjev, The yielding transition in periodically sheared binary glasses at finite temperature, Comput. Mater. Sci. 150, 162 (2018).
  * (28) A. D. S. Parmar, S. Kumar, and S. Sastry, Strain localization above the yielding point in cyclically deformed glasses, Phys. Rev. X 9, 021018 (2019).
  * (29) N. V. Priezjev, Slow relaxation dynamics in binary glasses during stress-controlled, tension-compression cyclic loading, Comput. Mater. Sci. 153, 235 (2018).
  * (30) N. V. Priezjev, Accelerated relaxation in disordered solids under cyclic loading with alternating shear orientation, J. Non-Cryst. Solids 525, 119683 (2019).
  * (31) N. V. Priezjev, Shear band formation in amorphous materials under oscillatory shear deformation, Metals 10, 300 (2020).
  * (32) W. T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki, and L. Berthier, Glass stability changes the nature of yielding under oscillatory shear, Phys. Rev. Lett. 124, 225502 (2020).
  * (33) N. V. Priezjev, Alternating shear orientation during cyclic loading facilitates yielding in amorphous materials, J. Mater. Eng. Perform. 29, 7328 (2020).
  * (34) N. V. Priezjev, A delayed yielding transition in mechanically annealed binary glasses at finite temperature, J. Non-Cryst. Solids 548, 120324 (2020).
  * (35) H. Bhaumik, G. Foffi, and S. Sastry, The role of annealing in determining the yielding behavior of glasses under cyclic shear deformation, PNAS 118, 2100227118 (2021).
  * (36) N. V. Priezjev, Accessing a broader range of energy states in metallic glasses by variable-amplitude oscillatory shear, J. Non-Cryst. Solids 560, 120746 (2021).
  * (37) N. V. Priezjev, Yielding transition in stable glasses periodically deformed at finite temperature, Comput. Mater. Sci. 200, 110831 (2021).
  * (38) S. Cui, H. Liu, and H. Peng, Anisotropic correlations of plasticity on the yielding of metallic glasses, Phys. Rev. E 106, 014607 (2022).
  * (39) P. Das, A. D. S. Parmar, and S. Sastry, Annealing glasses by cyclic shear deformation, J. Chem. Phys. 157, 044501 (2022).
  * (40) N. V. Priezjev, Mechanical annealing and yielding transition in cyclically sheared binary glasses, J. Non-Cryst. Solids 590, 121697 (2022).
  * (41) B. Shang, N. Jakse, P. Guan, W. Wang, and J.-L. Barrat, Influence of oscillatory shear on nucleation in metallic glasses: A molecular dynamics study, Acta Mater. 246, 118668 (2023).
  * (42) N. V. Priezjev, Fatigue failure of amorphous alloys under cyclic shear deformation, Comput. Mater. Sci. 226, 112230 (2023).
  * (43) V. V. Krishnan, K. Ramola, and S. Karmakar, Annealing effects of multidirectional oscillatory shear in model glass formers, Phys. Rev. Applied 19, 024004 (2023).
  * (44) N. V. Priezjev, Fatigue behavior of Cu-Zr metallic glasses under cyclic loading, Metals 13, 1606 (2023).
  * (45) L. Zella, J. Moon, and T. Egami, Ripples in the bottom of the potential energy landscape of metallic glass, Nat. Commun. 15, 1358 (2024).
  * (46) X. C. Tang and X. H. Yao, The effect of forced vibration coupling on amorphous alloy superplasticity, J. Non-Cryst. Solids 630, 122905 (2024).
  * (47) Y. Q. Cheng, E. Ma, and H. W. Sheng, Atomic level structure in multicomponent bulk metallic glass, Phys. Rev. Lett. 102, 245501 (2009).
  * (48) Y. Q. Cheng and E. Ma, Atomic-level structure and structure-property relationship in metallic glasses, Prog. Mater. Sci. 56, 379 (2011).
  * (49) Z. Fan and E. Ma, Predicting orientation-dependent plastic susceptibility from static structure in amorphous solids via deep learning, Nat. Commun. 12, 1506 (2021).
  * (50) S. J. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comp. Phys. 117, 1 (1995).
  * (51) K. Hyun, M. Wilhelm, C. O. Klein, K. S. Cho, J. G. Nam, K. H. Ahn, S. J. Lee, R. H. Ewoldt, and G. H. McKinley, A review of nonlinear oscillatory shear tests: Analysis and application of large amplitude oscillatory shear (LAOS), Prog. Polym. Sci. 36, 1697 (2011).
  * (52) M. L. Falk and J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids, Phys. Rev. E 57, 7192 (1998).
  * (53) X. Wang, W.-S. Xu, H. Zhang, and J. F. Douglas, Universal nature of dynamic heterogeneity in glass-forming liquids: A comparative study of metallic and polymeric glass-forming liquids, J. Chem. Phys. 151, 184503 (2019).