# Dynamic magnetostriction for antiferromagnets

Thomas Nussle<EMAIL_ADDRESS>CEA DAM/Le Ripault, BP 16, F-37260, Monts,
FRANCE  Institut Denis Poisson, Université de Tours, Université d’Orléans,
CNRS (UMR7013), Parc de Grandmont, F-37200, Tours, FRANCE  Pascal Thibaudeau
<EMAIL_ADDRESS>CEA DAM/Le Ripault, BP 16, F-37260, Monts, FRANCE
Stam Nicolis<EMAIL_ADDRESS>Institut Denis Poisson,
Université de Tours, Université d’Orléans, CNRS (UMR7013), Parc de Grandmont,
F-37200, Tours, FRANCE

###### Abstract

In this paper we study the switching properties of the dynamics of magnetic
moments, that interact with an elastic medium. To do so we construct a
Hamiltonian framework, that can take into account the dynamics in phase space
of the variables that describe the magnetic moments in a consistent way. It is
convenient to describe the magnetic moments as bilinears of anticommuting
variables that are their own conjugates. However, we show how it is possible
to avoid having to deal directly with the anticommuting variables themselves,
only using them to deduce non-trivial constraints on the magnetoelastic
couplings. We construct the appropriate Poisson bracket and a geometric
integration scheme, that is symplectic in the extended phase space and that
allows us to study the switching properties of the magnetization, that are
relevant for applications, for the case of a toy model for antiferromagnetic
NiO, under external stresses.

In the absence of magnetoelastic coupling, we recover the results reported in
the literature and in our previous studies. In the presence of the
magnetoelastic coupling, the characteristic oscillations of the mechanical
system have repercussions on the Néel order parameter dynamics. This is
particularly striking for the spin accumulation which is more than doubled by
the coupling to the strain ; here as well, the mechanical oscillations are
reflected on the magnetic dynamics. As a consequence of such a stress induced
strain, the switching time of the magnetization is slightly faster and the
amplitude of the magnetization enhanced.

##  I Introduction

With the emergence of spintronics, attention has now been focused on the
theoretical understanding and experimental manipulation of the properties of
magnetic materials at ever shorter length and time scales. In this context,
antiferromagnets (AFs) are the magnetically ordered materials which appear to
be among the most promising candidates for realizing fast spintronic devices
with very low power consumption Radu  _et al._ (2011);
*jungwirthAntiferromagneticSpintronics2016. Therefore quantitative modeling of
spin transport is also of topical interest.

Recent progress in describing spin transport and spin-transfert torque (STT)
effects in such AF devices opens a route towards developing multilevel memory
devices with switching speeds, that could exceed those of devices made of
ferromagnetic materials and semiconductors Železný  _et al._ (2018). Moreover
reversible strain-induced magnetization switching in ferromagnetic materials
has been reported that also allows the design of a rewritable, non-volatile,
non-toggle and extremely low energy straintronic memory Ahmad  _et al._
(2015). Recently a piezolectric, strain-controlled AF memory, insensitive to
magnetic fields, has been tested, as an example of controlling magnetism by
electric fields in multiferroic heterostructures Yan  _et al._ (2019).

However a theoretical description of all these phenomena is, still, very much,
work in progress, since, at the scales that are relevant, it is not possible
to separate the magnetic, electric and elastic responses–all must be treated
simultaneously.

It is, indeed, noteworthy that no direct local coupling between
electromagnetic field components is allowed in vacuum, at the classical level;
the magnetic response to an electric field is necessarily mediated by atomic
position arrangements of magnetic moments, which in response produce non-local
mechanical strains Fiebig (2005); Atulasimha and Bandyopadhyay (2010). This
response is known, more generally, under the heading of “magneto–elasticity”
Du Trémolet de Lacheisserie (1993).

Magneto-elastic phenomena have been, typically, viewed either from the
perspective of continuum mechanics, where the magnetic properties of materials
are incorporated into constitutive nonlinear laws of electro-conductive bodies
Maugin (1988), or through the introduction of effective magnetic anisotropies
mimicking the static strains of a solid Kittel (1949).

Most studies, however, do not investigate the backreaction of the
magnetization on the mechanical strains; however it is such a backreaction
that can become significant at very short length and time scales (where a
continuum description, anyway reaches its limits). This highlights the need
for developing multiscale models.

Even if multiscale models have been investigated for magneto-elastic couplings
Buiron  _et al._ (1999), with good experimental comparison for single crystals
and polycrystalline samples Daniel and Galopin (2008), insights from the
framework of localization and homogenization are still the mainstays for
deducing the constitutive laws for polycrystalline media. Moreover, in this
context, homogenization methods for multiscale mechanics assume the existence
of well–separated scales and different relations among length scales, that
lead to different effective equations, which, in turn, represent the
corresponding different physical effects, appropriate at each scale Mei and
Vernescu (2010).

Renormalization group ideas, of course, are instrumental in providing a
description of how the dynamics depends on the scale; what has been lacking,
to date, is a unified construction of the dynamics of elastic and spin degrees
of freedom, within a common Hamiltonian approach, at any particular scale. In
particular, what has, curiously, not been fully exploited is the resolution of
the constraints, that the spin degrees of freedom are subject to, since they
are their own canonically conjugate variables, as their equations of motion
are first order. This means that they can be most usefully expressed as
multilinear combinations of anticommuting variables. As we shall see, this can
lead to useful insights, in practice, about the dynamics of the coupled
system.

While the generalization of the canonical formalism, that incorporates the
constraints implied by spin degrees of freedom has been worked out a long time
ago Casalbuoni (1976a); *casalbuoniClassicalMechanicsBosefermi1976;
*berezinParticleSpinDynamics1977; *nelsonLagrangianTreatmentMagnetic1994, it
has not been applied to concrete situations of physical interest, simply
because these, up to now, were not of practical relevance.

In the present paper we provide a complete description of the combined
Hamiltonian dynamics of the elastic and of the spin degrees of freedom in
phase space and show how the anticommuting nature of the spin degrees of
freedom can be–indirectly–tested. It is the refined control over these and the
dynamical interplay between strain and magnetization that are the particular
novelties of our approach.

The elastic (local mechanical strain) and spin degrees of freedom are thus
assigned to sites of a lattice. We deduce the corresponding fully coupled
equations of motion from the corresponding Hamiltonian formalism, that treats
mechanical and spin degrees of freedom in a unified way and study the dynamics
of this system, comparing it to previous studies, first for the purely
magnetic part Cheng  _et al._ (2015) and then, also, for the magneto-
elastically coupled system Nussle  _et al._ (2019).

Of course to render the equations tractable we must resort to approximations.
The approximation we employ is the mean field approximation, that amounts, in
the present case, to a coarse–graining of the lattice, which means that the
sites are, in fact, domains, within each of which the properties are
homogeneous. Within this framework, our Hamiltonian formalism is, however,
exact.

The paper is organized as follows: In section II we construct a classical
Hamiltonian for the mechanical, and the spin degrees of freedom and discuss
their couplings. We discuss the advantages of describing the magnetic moments
in terms of anticommuting, Majorana, variables, whose salient properties are
reviewed in appendix A.

In section III, we introduce a Poisson bracket, which allows us to obtain the
equations of motion for the time-evolution of all the dynamical variables.
These are, then solved, in section IV using a symplectic integration scheme.
Special attention is devoted to the stability analysis of the schemes and to
the order of the different symplectic operators. This analysis is presented in
appendix B.

Finally in section V we illustrate the formalism for the case of a simple toy
model for NiO, consisting of two spins, interacting through antiferromagnetic
exchange, and switched by an external STT, as well as under an external
stress. It is here, in particular, that we show that the formalism is capable
of handling the interaction between two domains, taking into account the
coupling between elastic and magnetic degrees of freedom.

##  II The Hamiltonian

In this section we propose a Hamiltonian for the combined system of elastic
and spin degrees of freedom. These parts have, of course, been considered
before Du Trémolet de Lacheisserie (1993), so we shall review the salient
features, before introducing our approach.

The starting point for the elastic degrees of freedom is the framework of
mechanics for elastic solids, within the regime of the approximation of small
deformations, supplemented by the assumption of perfect mechanical micro-
reversibility Casimir (1945); *landauTheoryElasticity2008. These statements
imply that an elastic material is characterized by a collection of $N$ time-
dependant elastic deformations, which are described by a set of symmetric
Cauchy strain tensors ${\epsilon_{IJ}^{i}(t)}$ where $I,J$ are spatial
indices, ranging from 1 to 3, and assigned to each lattice site,
$i\in\\{1..N\\}$, of the material. At equilibrium, each of these variables is
independent of time and we assume that such a state not only exists, but can
be relevant for the time scales of interest Maugin (1988).

The total internal energy of this system can be described in terms of the
interactions between these variables, that take into account the elastic
mechanical part in the most general form and the interaction with external
stresses, which will be described by the magnetic degrees of freedom.
Schematically,

$\frac{\mathcal{H}}{V_{0}}=\frac{\mathcal{H}_{\mathrm{mech}}}{V_{0}}+\frac{\mathcal{H}_{\mathrm{ext}}}{V_{0}}$
(1)

We shall now spell out each term.

To lowest, non–trivial, i.e. quadratic, order,
$\mathcal{H}_{\mathrm{ext}}/V_{0}$ can be written as Gairola (1978);
*eringenNonlocalContinuumField2004

$\frac{{\cal{H}}_{\textrm{mech}}}{V_{0}}=\frac{1}{2}\sum_{i=1}^{N}C_{IJKL}^{i}\epsilon_{IJ}^{i}(t)\epsilon_{KL}^{i}(t),$
(2)

where an implicit sum on the repeated space indices is understood and where
$C_{IJKL}^{i}$ are the elastic constants with $V_{0}$ as the reference volume.

In this expression we have assumed that the different sites of the material do
not interact through elastic deformations (they will only interact through
spin exchange forces, to be spelled out below).

For an homogeneous material (i.e., $C^{(i)}=C$ for every site $i$), these
elastic constants can be expressed in terms of the two Lamé parameters Landau
_et al._ (2008), designated by $(C_{0},C_{1})$, which in Cartesian coordinates
read as follows:

$C_{IJKL}=C_{0}\delta_{IJ}\delta_{KL}+C_{1}\left(\delta_{IK}\delta_{JL}+\delta_{IL}\delta_{JK}\right).$
(3)

This mechanical system can be excited by an external stress
$\sigma_{IJ}^{\textrm{ext}}$, through the coupling

$\frac{{\cal
H}_{\textrm{ext}}}{V_{0}}=-\sum_{i=1}^{N}\sigma_{IJ}^{\textrm{ext}}\epsilon_{JI}^{i}.$
(4)

While this stress can be, in general, a time- and space-dependent quantity, it
is, nonetheless, assumed to be uniform over the sites of the material.

The magnetic degrees of freedom, noted $S_{I}^{i}(t)$, are assigned to the
same lattice sites as the strain and describe how their dynamics affects the
mechanical degrees of freedom, through the so-called magnetoelastic coupling.

The contribution of the magnetoelastic coupling to the Hamiltonian can be
deduced by assuming (global) Galilean invariance on the one hand and imposing
invariance under the symmetries of the point group Toupin and Bernstein
(1961); *callenStaticMagnetoelasticCoupling1963;
*brownMagnetoelasticInteractions1966 on the other. For an isotropic medium,
these requirements lead to the expressions for the total energy and the total
magnetization Rosen (1972). The total internal energy is, thus, a sum over the
lattice Tiersten (1964); *akhiezerSpinWaves1968.

Consequently, in the approximation of small deformations, the expansion of the
magnetoelastic coupling energy will contain even powers of the magnetization
${\bm{S}}^{i}$. It will, however, contain all powers of the strain tensor,
$\epsilon$.

The linear part in $\epsilon_{IJ}^{i}$ defines the first-order, magnetoelastic
coupling, which is responsible for the magnetostriction, while the quadratic
terms in $\epsilon_{IJ}^{i}$ define the second-order magnetoelastic coupling,
that is responsible for the “nonlinear response” in the elastic properties of
magnetic media, that have been, also, identified as describing “morphic
effects” Du Trémolet de Lacheisserie (1993).

Higher order terms are, usually, ignored, since their coefficients are assumed
to be much smaller than the lower order terms, already discussed.

Therefore, the contribution to the–total–Hamiltonian of the magnetoelastic
coupling can be written as

$\displaystyle{\cal H}_{\textrm{ME}}$
$\displaystyle=\sum_{i,j=1}^{N}B_{IJKL}^{(1)i,j}\epsilon_{IJ}^{i}(t)S_{K}^{i}(t)S_{L}^{j}(t)$
$\displaystyle+B_{IJKLMN}^{(2)i,j}\epsilon_{IJ}^{i}(t)\epsilon_{KL}^{j}(t)S_{M}^{i}(t)S_{N}^{j}(t)+\dots$
(5)

The total Hamiltonian can, thus, be written as

$\mathcal{H}_{\mathrm{tot}}=\mathcal{H}_{\mathrm{mech}}+\mathcal{H}_{\mathrm{ext}}+\mathcal{H}_{\mathrm{ME}}$
(6)

and leads to a “renormalization” of the effective elastic constants in the
following way: If one compares the last term of equation (II) with equation
(2), it is easy to see that a pair of effective elastic constants can be
introduced that depends on the value of the spins on each site, i.e.

$C_{IJLK}^{\textrm{eff},i,j}\equiv
C_{IJLK}^{ij}+\frac{1}{V_{0}}B_{IJKLMN}^{(2)i,j}S_{M}^{i}(t)S_{N}^{j}(t).$ (7)

This has been used to provide an interpretation of the “anomalous” temperature
dependence in the elastic constants in iron single crystals, for instance, as
resulting from the competition between spin ordering and diffusion effects
Dever (1972). In the present work however, we will only consider
$B_{IJKL}^{(1)i,j}$ and set aside the non-linear mechanical terms. Although
there is no particular difficulty to take the latter into account in this
formalism, we will take $B^{(2)}\equiv 0$ for the rest of this study and for
the sake of simplicity.

In Appendix A we recall that the most effective “classically equivalent
description” of the spin degrees of freedom, that can capture consistently
their interaction with the elastic degrees of freedom is not through the
variables ${\bm{S}}$ but through their “Doppelgänger” ${\bm{\xi}}$, related to
the $S_{I}^{i}$ through

$S_{I}^{i}\equiv-\frac{\imath}{2}\epsilon_{IJK}\xi_{J}^{i}\xi_{K}^{i}.$ (8)

The $\bm{\xi}$ can be identified with Majorana fermions, which have found many
applications recently in condensed matter systems Wilczek (2009);
*aliceaNewDirectionsPursuit2012, where new methods for controlling spin
degrees of freedom have been developed.

It is interesting to remark that the anticommuting variables, $\xi_{k}^{i}$
are not Grassmann variables, satisfying $\\{\xi_{I}^{i},\xi_{J}^{j}\\}=0$;
but, rather, $\\{\xi_{I}^{i},\xi_{J}^{i}\\}=\delta^{ij}\delta_{IJ}$, i.e. that
the $\xi_{I}$ generate a Clifford algebra on each site Lounesto (2003). It is
in this way that the $S_{I}^{i}$, defined through eq. (8), satisfy the angular
momentum algebra Varshalovich  _et al._ (1988).

If one were tempted to simply replace ${\bm{S}}$ by ${\bm{\xi}}$ in eq.(II),
it is interesting to remark that, for $N=1$, and because of the symmetries of
$B^{(1)}$ as recalled in reference Du Trémolet de Lacheisserie (1993), ${\cal
H}_{ME}=0$, which implies that the dynamics as it is cannot be encoded by this
Hamiltonian.

On the other hand, this allows us to understand the constraints on the allowed
terms in the true Hamiltonian. They must, necessarily, involve more than one
spins. Indeed, we can construct expressions that are multi-linear combinations
of the $S^{i}_{I}$ on different sites, potentially, up to order $N$, the
number of sites, since no two identical $\xi_{I}^{i}$ variables occur in the
same monomial. This is, therefore, a nice way of automatically organizing the
multi-spin terms of the Hamiltonian.

In the formalism of Atomistic Spin Dynamics, since the magnetic moments are
localized Evans  _et al._ (2014); Eriksson  _et al._ (2017), it is customary
to consider spatial averages, around each site, defining an effective
macroscopic localized spin

$S_{I}^{\textrm{eff}}=\langle S_{I}^{i}\rangle_{i}$ (9)

which implies that multi-linear expressions, $S_{I}S_{J}$ no longer vanish
identically, as was imposed by the anticommuting nature of $\bm{\xi}$ before.
We can, thus, understand the relevance of this averaging procedure, in terms
of the description of the spin degrees of freedom in terms of the
anticommuting variables.

The justification for equation (II) has been advocated a long time ago van
Vleck (1937); *neelAnisotropieMagnetiqueSuperficielle1954;
*leeMagnetostrictionMagnetomechanicalEffects1955, as stemming from a model for
a two–body interaction, that is itself a pedagogical version of the quantum
theory of interacting magnons and phonons Sabiryanov and Jaswal (1999). In eq.
(II) what has been left unspecified are the properties of the tensor $B^{(1)}$
under exchange of the indices $i$ and $j$, that label the sites. What is
customary is to assume that $B^{(1)}$, in fact, does not depend on $i$ and $j$
at all–it is homogeneous across the material. By stopping the expansion at
first order in the mechanical deformation in (II) and assuming homogeneity in
the material, we, therefore, end up with the following expression

${\cal H}_{\textrm{ME}}=B_{IJKL}^{(1)}\epsilon_{IJ}(t)S_{K}(t)S_{L}(t).$ (10)

Again, in the case of an isotropic medium, the elements in $B^{(1)}$ enjoy the
same symmetries than the elastic constants and can be expressed only in terms
of two constants $B^{(1)}_{0}$ and $B^{(1)}_{1}$ by the following expression

$B_{IJKL}^{(1)}=B^{(1)}_{0}\delta_{IJ}\delta_{KL}+B^{(1)}_{1}\left(\delta_{IK}\delta_{JL}+\delta_{IL}\delta_{JK}\right).$
(11)

In summary, the microscopic theory underlying Eq.(10) describes the
interaction of three particles, two of which describe spin states, and one the
state of the elastic deformation. However the spin states, are, in fact, bound
states of more “fundamental” entities, the anticommuting $\bm{\xi}$. $B^{(1)}$
can then be interpreted as the vertex of this interaction, that describes a
spinning particle, that does not directly interact with itself–since any such
interaction is inconsistent with the fact that the $\bm{\xi}$ anticommute.
Such a self-interaction can, however, appear on larger scales, when spatial
averages can become meaningful for describing the spin degrees of freedom.

In conclusion to this section, the Hamiltonian we have constructed describes
the interaction of magnetic moments through their embedding within an elastic
medium. What remains to be done is to define the Poisson brackets, that can
take into account the evolution in phase space of commuting, as well as
anticommuting degrees of freedom. Indeed, the construction procedure for
commuting canonical variables is well known, however, including anticommuting
variables in a unified way has remained a rather esoteric subject–known in
theory Casalbuoni (1976a, b); Berezin and Marinov (1977), not, however,
implemented, in practice.

In the following section, we shall construct the equations of motion in the
phase space of the elastic and the spinning degrees of freedom, that
implements these ideas in practice.

In order to work directly with the anticommuting variables themselves, in
combination with the mechanical degrees of freedom requires implementing a
“graded Poisson bracket”; one way to do this is discussed in Appendix A.

##  III The Poisson brackets and the equations of motion

To deduce from the Hamiltonian, discussed in the previous section, the
equations of motion, we must define the appropriate pairs of canonically
conjugate variables and consequently their Poisson bracket.

First we recognize that $\varepsilon$ acts as a tensor and one can understand
the definition of the Poisson bracket of rank-2 symmetric tensors as an
application of the DeDonder-Weyl covariant hamiltonian formulation of field
theory Kastrup (1983); *kanatchikovCanonicalStructureDonderWeyl1993. Although
the context is different, the ADM procedure in general relativity also
provides such a Poisson bracket Arnowitt  _et al._ (1960) (with further
relations, between the conjugate variables that are not relevant here).

Although no clear consensus has, in fact, emerged on the properties of Poisson
bracket of rank-2 tensors Hojman  _et al._ (1976), if one focuses on the
special case of strain tensors, that depend only on time, the following
conjugate momentum can be introduced

$\pi_{IJ}\equiv\frac{\partial\mathscr{L}}{\partial{\dot{\varepsilon}_{IJ}}},$
(12)

where $\mathscr{L}(\varepsilon_{IJ}(t),\dot{\varepsilon}_{IJ}(t))$ is the
unconstrained and free Lagrangian. $\dot{\varepsilon}_{IJ}(t)$ are the
components of the strain-rate tensor Hill (1970).

Thus, we can build the corresponding Hamiltonian ${\cal H}$ for the time
evolution with tensor variables for mechanical deformations $\bm{\varepsilon}$
and their conjugated momenta $\bm{\pi}$, because these quantities admit
unbounded numerical values, as the corresponding Legendre transform

${\cal
H}({\bm{\varepsilon}},{\bm{\pi}})=\pi_{IJ}\dot{\varepsilon}_{IJ}-\mathscr{L},$
(13)

up to a total time derivative for ${\mathscr{L}}$.

For the mechanical system only (i.e. for functions
$A({\bm{\varepsilon}},{\bm{\pi}})$ and $B({\bm{\varepsilon}},{\bm{\pi}})$),
the Poisson bracket can be defined in perfect analogy to that of any particle
system, that explores a given target space geometry (to which refer the
indices $I,J,K,L$), by the usual relations

$\\{A,B\\}_{PB}=\frac{\partial A}{\partial\varepsilon_{IJ}}\frac{\partial
B}{\partial\pi_{IJ}}-\frac{\partial A}{\partial\pi_{IJ}}\frac{\partial
B}{\partial\varepsilon_{IJ}}.$ (14)

In our case, the dynamical variables are the real symmetric rank-two tensors,
$\varepsilon_{IJ}$ and $\pi_{IJ}$ ($I,J=1,2,3$), which are canonically
conjugate in the sense that their Poisson brackets are deemed to satisfy the
following properties:

$\displaystyle\\{\varepsilon_{IJ},\pi_{KL}\\}_{PB}$
$\displaystyle=\delta_{IJKL},$ (15)
$\displaystyle\\{\varepsilon_{IJ},\varepsilon_{KL}\\}_{PB}$
$\displaystyle=\\{\pi_{IJ},\pi_{KL}\\}_{PB}=0$ (16)

where $\delta_{IJKL}$ is a $\delta$ “tensor”, reflecting the real, symmetric,
nature of the Poisson brackets Kijowski and Szczyrba (1976);
*szczyrbaSymplecticStructureSet1976; *marsdenCovariantPoissonBrackets1986 and
defined as a product of Kronecker $\delta$s. Very schematically, we may write

$\delta_{IJKL}=\left\\{\begin{array}[]{rcl}\delta_{IJ}\delta_{KL}\\\
\delta_{IJ}\delta_{LK}\\\ \delta_{JI}\delta_{LK}\\\ \delta_{JI}\delta_{KL}\\\
\end{array}\right.$ (17)

where each choice of the RHS corresponds to a choice of indices in the Poisson
brackets. This choice can be supplemented by any linear combination of these
$\delta$s that enforces the symmetries of the tensors.

We now wish to include as phase space coordinates, the components of the spin
vector, $\bm{S}$. We follow reference Yang and Hirschfelder (1980) and the
details are summarized in appendix A. The “generalized” Poisson bracket, for
the canonical variables of our system, can be written as

$\begin{array}[]{rcl}\displaystyle\\{A,B\\}_{PB}\equiv\frac{\partial
A}{\partial\varepsilon_{IJ}}\frac{\partial
B}{\partial\pi_{IJ}}&-&\displaystyle\frac{\partial
A}{\partial\pi_{IJ}}\frac{\partial B}{\partial\varepsilon_{IJ}}\\\
&-&\displaystyle\frac{1}{\hbar}\epsilon_{IJK}S_{I}\frac{\partial A}{\partial
S_{J}}\frac{\partial B}{\partial S_{K}},\end{array}$ (18)

where $\hbar$ is introduced to restore the physical dimensions of the Poisson
bracket, since ${\bm{S}}$ is considered dimensionless, for the sake of
simplicity.

It should be stressed that this $\hbar$ does not imply that any quantum
effects are present, since the dynamics is fully classical. It is simply a
bookkeeping device for a quantity that has the dimensions of angular
momentum–i.e. of an area in phase space–and reflects the fact that the
equation of motion for $S_{I}$ is of first order. Quantum fluctuations will
introduce the “real” $\hbar$.

Using this Poisson bracket, we can obtain the equations of motion for the
phase space variables:

$\begin{array}[]{rcl}\displaystyle\dot{\varepsilon}_{IJ}&=&\displaystyle\\{\varepsilon_{IJ},{\cal
H}\\}_{PB}=\frac{\partial{\cal H}}{\partial\pi_{IJ}}\\\
\displaystyle\dot{\pi}_{IJ}&=&\displaystyle\\{\pi_{IJ},{\cal
H}\\}_{PB}=-\frac{\partial{\cal H}}{\partial\varepsilon_{IJ}}\\\
\displaystyle\dot{S}_{I}&=&\\{S_{I},{\cal
H}\\}_{PB}=\displaystyle\frac{1}{\hbar}\varepsilon_{IJK}S_{J}\frac{\partial{\cal
H}}{\partial S_{K}}\end{array}$ (19)

The consistency of this formalism can be checked by noting that these
equations preserve the volume in phase space

$\frac{\partial\dot{\varepsilon}_{IJ}}{\partial\varepsilon_{IJ}}+\frac{\partial\dot{\pi}_{IJ}}{\partial\pi_{IJ}}+\frac{\partial\dot{S}_{K}}{\partial
S_{K}}=0.$ (20)

That the dynamics preserves the volume in phase space does not, of course,
imply anything about whether the system thus defined is integrable or shows
Hamiltonian chaos.

The internal energy $U$, involving the mechanical energy for the deformed
elastic medium Landau  _et al._ (2008), the magnetic energy, defined by the
Zeeman term Skomski (2008) and the magneto-elastic energy Du Trémolet de
Lacheisserie (1993), that takes into account the interaction of the magnetic
moment with the medium, takes the form

$\displaystyle U=$
$\displaystyle\frac{V_{0}}{2}C_{IJKL}\varepsilon_{IJ}\varepsilon_{KL}-V_{0}\sigma^{\textrm{ext}}_{IJ}\varepsilon_{IJ}+B_{IJKL}^{(1)}\varepsilon_{IJ}S_{K}S_{L}$
$\displaystyle-\hbar\omega_{I}S_{I},$ (21)

where $C$ is the fully symmetric tensor, defining the elastic response,
$\sigma^{\textrm{ext}}$ is the external stress tensor, ${\bm{\omega}}$ is the
effective external magnetic field (expressed as a frequency) and $B^{(1)}$ is
the fully symmetric linear magnetostriction tensor.

The “kinetic” term, containing the conjugate momenta, can be written,
schematically, as

${\cal H}_{\mathrm{kinetic}}=\frac{1}{2}\pi_{IJ}M^{-1}_{IJKL}\pi_{KL}$ (22)

where $M$ is a fully symmetric “mass” matrix, that describes the inertial
response.

For the case of isotropic materials, it is assuming that the $M$ tensor has
the form given by Lamé, with only two characteristic constants:

$M_{IJKL}=M_{0}\delta_{IJ}\delta_{KL}+M_{1}\left(\delta_{IK}\delta_{JL}+\delta_{IL}\delta_{JK}\right).$
(23)

The tensors $C$ and $B$ are decomposed in the same way.

Consequently, the inverse of these tensors can then be deduced from
$M_{IJKL}M^{-1}_{IJMN}=\frac{1}{2}\left(\delta_{KM}\delta_{LN}+\delta_{KN}\delta_{LM}\right)$.
Thus

$\begin{array}[]{rcl}M^{-1}_{IJKL}&=&\displaystyle\frac{-M_{0}}{2M_{1}(3M_{0}+2M_{1})}\delta_{IJ}\delta_{KL}\\\
&+&\displaystyle\frac{1}{4M_{1}}\left(\delta_{IK}\delta_{JL}+\delta_{IL}\delta_{JK}\right)\end{array}$
(24)

and the equations of motion become

$\left\\{\begin{array}[]{rcl}\dot{\epsilon}_{IJ}&=&M^{-1}_{IJKL}\pi_{KL}\\\
\dot{\pi}_{IJ}&=&-V_{0}C_{IJKL}\epsilon_{KL}+V_{0}\sigma^{\textrm{ext}}_{IJ}-B_{IJKL}S_{K}S_{L}\\\
\dot{S}_{I}&=&\varepsilon_{IJK}\left(\omega_{J}-\frac{2}{\hbar}B_{ABJC}^{(1)}\epsilon_{AB}S_{C}\right)S_{K}\end{array}\right.$
(25)

highlighting how the mechanical and magnetic subsystems are coupled.

The last equation–as expected!–is a precession equation for the components of
${\bm{S}}$ around both the effective field ${\bm{\omega}}$ and an additional
field, that depends on the strain tensor and the spin vector.

In the following section we shall show how to solve these equations, in a way
that preserves the symmetries of the extended phase space.

##  IV Geometric integration

Solving the coupled system of equations (25) is the next step.

Since, in the previous section, we have shown that these equations describe a
volume preserving transformation of the enlarged phase space, encompassing
elastic and spin variables, it is natural to rewrite them, in terms of the
action of a Liouville operator. Therefore, we shall write eqs.(25) as

$\begin{array}[]{rcl}{\dot{\bm{\varepsilon}}}&=&{{\cal
L}}_{\bm{\epsilon}}{\bm{\varepsilon}},\\\ {\dot{\bm{\pi}}}&=&{\cal
L}_{\bm{\pi}}{\bm{\pi}},\\\ {\dot{\bm{S}}}&=&{\cal
L}_{\bm{S}}{\bm{S}}.\end{array}$ (26)

where ${\cal L}$ is the Liouville operator. This formulation allows us to
implement, manifestly, time-reversible, area preserving algorithms, for
solving these equations numerically.

The general scheme is as follows: Consider an arbitrary function $f$ of the
canonically conjugate variables of our many-body system. This function,
$f({\bm{\varepsilon}},{\bm{\pi}},{\bm{S}})$, depends on the time $t$
implicitly; that is, through the dependence of
$({\bm{\varepsilon}},{\bm{\pi}},{\bm{S}})$ on $t$. The time derivative of $f$
is $\dot{f}$ such as

$\displaystyle{\dot{f}}$ $\displaystyle=\dot{\varepsilon}_{IJ}\frac{\partial
f}{\partial\varepsilon_{IJ}}+\dot{\pi}_{IJ}\frac{\partial
f}{\partial\pi_{IJ}}+\dot{S}_{I}\frac{\partial f}{\partial S_{I}}$
$\displaystyle\equiv{\cal L}f.$ (27)

The last line defines the (total) Liouville operator

${\cal
L}=\dot{\varepsilon}_{IJ}\frac{\partial}{\partial\varepsilon_{IJ}}+\dot{\pi}_{IJ}\frac{\partial}{\partial\pi_{IJ}}+\dot{S}_{I}\frac{\partial}{\partial
S_{I}}.$ (28)

Equation (27) can be integrated formally as an initial value problem to obtain
$f$ at any time :

$f({\bm{\varepsilon}}(t),{\bm{\pi}}(t),{\bm{S}}(t))=e^{{\cal
L}t}f({\bm{\varepsilon}}(0),{\bm{\pi}}(0),{\bm{S}}(0)).$ (29)

It is not difficult to see that ${\cal L}={\cal L}_{\bm{\varepsilon}}+{\cal
L}_{\bm{\pi}}+{\cal L}_{\bm{S}}$. However, these single Liouville operators do
not commute two-by-two, as the reader may easily check by computing ${\cal
L}_{\bm{u}}{\cal L}_{\bm{v}}f-{\cal L}_{\bm{v}}{\cal L}_{\bm{u}}f\neq 0$ for
any function $f$ and any combination $(u,v)$ of the individual Liouville
operator ${\cal L}_{\bm{\varepsilon}},{\cal L}_{\bm{\pi}},{\cal L}_{\bm{S}}$.
This means that

$e^{{\cal L}t}=e^{{\cal L}_{\bm{\varepsilon}}t+{\cal L}_{\bm{\pi}}t+{\cal
L}_{\bm{S}}t}\neq e^{{\cal L}_{\bm{\varepsilon}}t}e^{{\cal
L}_{\bm{\pi}}t}e^{{\cal L}_{\bm{S}}t}.$ (30)

According to the Magnus expansion Blanes  _et al._ (2009), however, it is
always possible to express $e^{{\cal L}t}$ as a product of the individual
operators at any given order in time, according to the so–called “splitting
method” Hairer  _et al._ (2006). This ensures that the numerical algorithm
preserves phase space volumes exactly.

For instance, for a fixed timestep $\tau$, upon expanding up to the third
order in time, the following sequence of products

$\displaystyle e^{{\cal L}\tau}$ $\displaystyle=e^{{\cal
L}_{\bm{S}}\frac{\tau}{4}}e^{{\cal L}_{\bm{\pi}}\frac{\tau}{2}}e^{{\cal
L}_{\bm{S}}\frac{\tau}{4}}e^{{\cal L}_{\bm{\varepsilon}}\tau}e^{{\cal
L}_{\bm{S}}\frac{\tau}{4}}e^{{\cal L}_{\bm{\pi}}\frac{\tau}{2}}e^{{\cal
L}_{\bm{S}}\frac{\tau}{4}}+{\cal O}(\tau^{3}),$ (31)

can be generated; this sequence is, in fact, one of six possible, that have
the property of preserving the symplectic structure of the Poisson brackets.
Therefore, any one of them can be chosen. The possible combinations are
presented in table 1.

$\begin{array}[]{c|ccccccc}A&{\frac{\tau}{4}}&{\frac{\tau}{2}}&{\frac{\tau}{4}}&{\tau}&{\frac{\tau}{4}}&{\frac{\tau}{2}}&{\frac{\tau}{4}}\\\
\hline\cr
1&{\bm{S}}&{\bm{\pi}}&{\bm{S}}&{\bm{\epsilon}}&{\bm{S}}&{\bm{\pi}}&{\bm{S}}\\\
2&{\bm{S}}&{\bm{\epsilon}}&{\bm{S}}&{\bm{\pi}}&{\bm{S}}&{\bm{\epsilon}}&{\bm{S}}\\\
3&{\bm{\pi}}&{\bm{\epsilon}}&{\bm{\pi}}&{\bm{S}}&{\bm{\pi}}&{\bm{\epsilon}}&{\bm{\pi}}\\\
4&{\bm{\pi}}&{\bm{S}}&{\bm{\pi}}&{\bm{\epsilon}}&{\bm{\pi}}&{\bm{S}}&{\bm{\pi}}\\\
5&{\bm{\epsilon}}&{\bm{\pi}}&{\bm{\epsilon}}&{\bm{S}}&{\bm{\epsilon}}&{\bm{\pi}}&{\bm{\epsilon}}\\\
6&{\bm{\epsilon}}&{\bm{S}}&{\bm{\epsilon}}&{\bm{\pi}}&{\bm{\epsilon}}&{\bm{S}}&{\bm{\epsilon}}\\\
\end{array}$ Table 1: Decomposition table of symplectic integrators

While these schemes are free from “global” errors, they are, of course,
sensitive to “local” errors, due to the finite value of the timestep. It is,
also, not at all obvious that all six can be implemented with comparable
efficiency. It is, therefore, useful to study the numerical stability and
efficiency of these different combinations, in particular, as former studies
in molecular dynamics Batcho and Schlick (2001) and magnetic molecular
dynamics Beaujouan  _et al._ (2012) showed, apparently, numerical differences
between them.

A sampler of such a study is presented in appendix B.

It is important to keep in mind that the one–step evolution operators for
${\bm{\varepsilon}}$ and ${\bm{\pi}}$ describe shifts of the corresponding
tensor components, whereas the one–step evolution operator for ${\bm{S}}$
describes rotations. In equations

$\displaystyle e^{{\cal
L}_{\varepsilon}\tau}({\bm{\varepsilon}}(t),{\bm{\pi}}(t),{\bm{S}}(t))$
$\displaystyle=({\bm{\varepsilon}}(t)+\tau{\dot{\bm{\varepsilon}}(t)},{\bm{\pi}}(t),{\bm{S}}(t))$
(32) $\displaystyle e^{{\cal
L}_{\pi}\tau}({\bm{\varepsilon}}(t),{\bm{\pi}}(t),{\bm{S}}(0))$
$\displaystyle=({\bm{\varepsilon}}(t),{\bm{\pi}}(t)+\tau{\dot{\bm{\pi}}(t)},{\bm{S}}(t))$
(33) $\displaystyle e^{{\cal
L}_{S}\tau}({\bm{\varepsilon}}(t),{\bm{\pi}}(t),{\bm{S}}(t))$
$\displaystyle=({\bm{\varepsilon}}(t),{\bm{\pi}}(t),\underbrace{{\bm{S}}(t+\tau)}_{R(\tau){\bm{S}}(t)})$
(34)

where ${\bm{S}}(t+\tau)=R(\tau){\bm{S}}(t)$ is given by the Rodrigues’
rotation formula Rodrigues (1840); *honerkampTheoreticalPhysicsClassical1993;
*thibaudeauThermostattingAtomicSpin2012 for a spin vector around a given
rotation vector ${\bm{\tilde{\omega}}}(t)$ where each of its components are
$\tilde{\omega}_{I}(t)=\omega_{I}(t)-\frac{2}{\hbar}B^{(1)}_{JKLI}\varepsilon_{JK}(t)S_{L}(t)$.
These equations describe the phase space of one particle only. To describe the
dynamics of a continuum, we must deduce the equations for many particles.

One way to generalize eqs.(25) for the case of many particles, labeled by an
index, $i=1,\ldots,N$, according to the conventions of the previous sections,
is the following

$\displaystyle\dot{\epsilon}^{i}_{IJ}$
$\displaystyle=[M^{i}]^{-1}_{IJKL}\pi^{i}_{KL}$ (35)
$\displaystyle\dot{\pi}^{i}_{IJ}$
$\displaystyle=-V_{0}C^{i}_{IJKL}\epsilon^{i}_{KL}+V_{0}\sigma^{i\hskip
2.0pt\textrm{ext}}_{IJ}-B^{i}_{IJKL}S^{i}_{K}S^{i}_{L}$ (36)
$\displaystyle\dot{S}^{i}_{I}$
$\displaystyle=\varepsilon_{IJK}\left(\omega^{i}_{J}-\frac{2}{\hbar}B^{(1)i}_{ABJC}\epsilon^{i}_{AB}S^{i}_{C}\right)S^{i}_{K}$
(37)

The case of a staggered AF state is treated simply by letting $N=2$. In order
to simplify the mechanical part further, we can impose additional conditions
pertaining to the uniformity of the external stress, mechanical constants,
mass matrices and magneto-elastic constants, at each site. In what follows, we
shall use the following Ansatz:

$\begin{array}[]{rcl}B^{(1)1}_{IJKL}&=&B^{(1)2}_{IJKL},\\\
C^{1}_{IJKL}&=&C^{2}_{IJKL},\\\ M^{1}_{IJKL}&=&M^{2}_{IJKL},\\\
\sigma^{1\hskip 2.0pt\textrm{ext}}_{IJ}&=&\sigma^{2\hskip
2.0pt\textrm{ext}}_{IJ}.\end{array}$ (38)

The conservative part of the precession contains a local field, which is
modified to include the antiferromagnetic exchange between the sites and a
single anisotropy axis ${\bm{n}}$ with an intensity $K$. One has

$\displaystyle\omega_{I}^{i}$
$\displaystyle=\frac{1}{\hbar}\sum_{<ij>}J^{ij}S_{I}^{j}+\frac{K}{\hbar}n_{J}s_{J}^{i}n_{I}$
(39)

Because of the exchange field, the Liouville operators for different spins do
not commute either. A global geometric integrator, implementing the approach
of Omelyan Omelyan  _et al._ (2003), which remains accurate up to third order
in the timestep expansion must, therefore, be constructed.

For any given timestep, $\tau$, the corresponding expression for the evolution
operator, reads

$\displaystyle e^{{\cal L}_{S}\tau}$ $\displaystyle=e^{{\cal
L}_{S^{1}}\frac{\tau}{2}}e^{{\cal L}_{S^{2}}\tau}e^{{\cal
L}_{S^{1}}\frac{\tau}{2}}+{\cal O}(\tau^{3}).$ (40)

For $N=2$, this operator is numerically identical to the operator, obtained by
permuting the site indices, $1\leftrightarrow 2$.

The same reasoning is applied for the Liouville operators for the elastic
variables, that enter in equations (26) and the corresponding global geometric
integrators are constructed along the same lines. The system is then
integrated by following one of the schemes displayed in table 1.

With these tools, we can study a plethora of phenomena, that are sensitive to
the coupling of magnetic, electric and elastic degrees of freedom. In the
following section we shall apply this formalism for studying “switching”
effects of the magnetization in antiferromagnetic media. The numerical
accuracy ensured by the geometric integrators is necessary to describe
picosecond switching times.

##  V Antiferromagnetic ultrafast switching under stress

In this section, we shall apply the formalism constructed above, to describe
how it is possible to generate and manipulate picosecond switching of the
magnetization, induced by STTs, from a short pulse of electric current, in
elastic media that exhibit staggered AF order.

Figure 1: (Color online) Switching scheme for antiferromagnetic magneto-
mechanically coupled toy model with external STT.

Such a fast switching process in AFs is schematically displayed in figure 1
and can be realized, in practice, by a femtosecond laser excitation of the
magnetic moments that generate a field-induced STT on the two sub-lattices.
During the pulse, energy is transferred from conduction electrons to the sub-
lattice magnetic moments via STT (this is how the electric current affects the
magnetic response), and hence contributes to the exchange energy between the
moments, since local moments can be canted noncollinearly. The energy due to
the strong exchange interaction between neighboring moments, which is commonly
found in AFs, in particular due to magnetic anisotropy, appears as an
effective inertia to the motion, leading to the appearance of a time scale,
much longer than that of the excitation pulse. Afterwards, the system follows
a natural path along the easy-plane to circumvent the unfavorable anisotropy
barrier and finally relaxes to a new magnetic configuration. The resulting
dynamics is a switching of the two sub-lattices moments to the opposite
direction through the easy-plane.

This process defines the ultrafast antiferromagnetic switching effect. During
the process, the Néel vector, that probes the difference between the magnetic
moments of the two sub-lattices, acquires a net value that can be transferred
as a so-called “spin accumulation”, to an adjacent non-magnetic normal
material in order to pump the produced spin current via scattering of
electrons. The reverse mechanism can be realized, as well.

What was not considered before is the possibility to enhance or inhibit such a
switching, depending on tensile or compression effects, that are generated by
an external stress, that couples to the internal strain, produced by the
intrinsic magneto-elastic interaction in such materials, as depicted in figure
1.

In order to conform to the notation used in our previous studies Nussle  _et
al._ (2019), equation (37) is supplemented with a non-conservative part,
labeled ${\bm{T}}$, on the RHS, which includes both a transverse damping and a
damping-like STT torque

$\displaystyle T_{I}^{i}$
$\displaystyle=\alpha\epsilon_{IJK}S_{J}^{i}\dot{S}_{K}^{i}+G\left(s_{I}^{i}s_{J}^{i}p_{J}-p_{I}s_{J}^{i}s_{J}^{i}\right).$
(41)

The resulting equations of motion, are numerically integrated, using the
approach of section V. It is noteworthy that, this equation, also, describes
the “backreaction” of the magnetic response on the spin transfer torque. On
the other hand, we do not consider how this backreaction affects the current
pulse itself, that is assumed external.

In figures 2 and 3 we report the evolution of the average magnetization
${\bm{m}}\equiv\frac{1}{2}\left({\bm{S}}^{1}+{\bm{S}}^{2}\right)$ and the Néel
vector ${\bm{l}}\equiv\frac{1}{2}\left({\bm{S}}^{1}-{\bm{S}}^{2}\right)$, in
the presence as well as the absence of magnetoelastic coupling, for a moderate
external stress of $\sigma_{xx}=-2\mu_{0}M_{s}^{2}$.

Figure 2: (Color online) Average magnetization (upper panel) and Néel order
parameter components (middle panel) for uncoupled switching. {$K=2\pi$rad.GHz,
$\omega_{E}$=172.16 rad.THz, $M_{s}=5.10^{5}\mathrm{A.m}^{-1}$}. Initial
conditions: {${\bm{s}}^{1}(0)=-{\bm{s}}^{2}(0)=\hat{\bm{x}}$}. The lower panel
displays the STT pulses. The figures agree with the referenceCheng  _et al._
(2015) because the magnetoelastic constants are set to zero. Figure 3: (Color
online) Average magnetization (upper panel) and Néel order parameter
components (middle panel) for a coupled switching. {$K=2\pi$rad.GHz,
$\omega_{E}$=172.16 rad.THz, $M_{s}=5.10^{5}\mathrm{A.m}^{-1}$}. Initial
conditions: {${\bm{s}}^{1}(0)=-{\bm{s}}^{2}(0)=\hat{\bm{x}}$,
$\epsilon^{i}_{IJ}(0)=0$}. The lower panel displays the STT pulses. Here
$B_{0}^{(1)}=7.7\mu_{0}M_{s}^{2}$ and $B_{1}^{(1)}=-23\mu_{0}M_{s}^{2}$.

We start the simulations using an initial configuration, where spins are
aligned along the $\hat{\bm{x}}$-axis in an antiferromagnetic configuration
and apply two $10$ ps electric pulses of $0.0034$ rad.THz intensity, each
separated by $50$ ps, in the ${\bm{p}}={\hat{\bm{z}}}$-direction. In addition
to the exchange interaction, the spins are subjected to a global anisotropy,
along the ${\bm{n}}={\hat{\bm{x}}}$-axis. The numerical value for $\alpha$ in
the following studies is $0.005$.

In the absence of magnetoelastic coupling, our results are identical to those
by Cheng et al. Cheng  _et al._ (2015) and to those we obtained under the same
conditions in our previous work Nussle  _et al._ (2019).

In the presence of magnetoelastic coupling, however, as the mechanical system
is undamped, the results are quite different. Because of the presence of a
constant stress, a finite mass matrix and non-zero elastic constants, the
mechanical system is expected to be oscillating freely, which is, indeed,
observed in figure 4, where we display the evolution of the strain components
over time. Indeed, as one can see on figures 3 and 4 the characteristic
oscillations of the mechanical system have repercussions on the Néel order
parameter dynamics. This is even more striking for the spin accumulation which
is more than doubled by the coupling to the strain ; here as well, the
mechanical oscillations are reflected on the magnetic dynamics. What is
interesting is that these characteristic mechanical oscillations are expected
to be related to the sound velocity of the medium Ogden (1997), which is
controlled by the constants $M_{0}$ and $M_{1}$. In the present study, we have
chosen $M_{0}=0$ and $M_{1}=10000$, but there should be a way to obtain these
coefficients of the mass matrix more efficiently, in order to describe these
oscillations more accurately from experimentally easily accessible data. This,
however is not the focus of the present study and shall be studied more
thoroughly elsewhere.

Figure 4: (Color online) Strain components as a function of time with the
magnetoelastic constants turned on with parameters identical to figure 3.

As a consequence of such a stress induced strain, the switching time of the
magnetization is slightly faster and the amplitude of the magnetization
enhanced, depending on the values of the stress and magnetoelastic constants,
as depicted in figure 3. We observe that a compression in one direction
becomes tensile in the other one and vice-versa depending on the sign of
$B_{1}^{(1)}$. Conversely, if $B_{1}^{(1)}$ is of the opposite sign, a tensile
stress would slow the switching process down.

This leads to another interesting question which is how the switching time is
affected by the values taken by either $B_{0}^{(1)}$ and $B_{1}^{(1)}$. We
have verified that changing $B_{0}^{(1)}$ has no influence on it because this
coefficient does not enter in the precession equation. Its simply produces the
Joule magnetostriction phenomena as anticipated before in reference Du
Trémolet de Lacheisserie (1993). However, changing $B_{1}^{(1)}$ has a strong
influence. The Figure 5 displays the switching time as a function of the
values of varying $B_{1}^{(1)}$. By defining $R$ as the ratio of $B_{1}^{(1)}$
divided by the its natural value in NiO, we perfom different simulations and
report the time where the Neel vector $x$-component crosses over to negative
values. One can see that this time decreases by increasing $B_{1}^{(1)}$ which
makes sense, as this is amplifying the effects of the tensile stress. When $R$
becomes lower than $1$, the effect is opposite.

Figure 5: Switching time (in ps) as a function of varying $B_{1}^{(1)}$ over
its natural value in NiO.

One aspect which is however surprising is that there is a quick drop of the
switching time for approximately $4$ times the value of $B_{1}^{(1)}$ over its
natural value in NiO. This is due to the oscillations of the magnetization
induced by the elastic coupling. As this coupling becomes stronger, one of the
mechanical oscillatory peaks dips below 0 and thus the switching time appears
to drop in a discontinuous fashion. This is another interesting consequence to
the interplay of the magnetization and the strain, that could be captured in
real experiments.

##  VI Discussion

In the present paper we have presented the details for a Hamiltonian framework
that can describe the dynamics, as well as the coupling, between magnetism and
elasticity by a careful definition of the corresponding dynamical variables in
an extended phase-space.

The magnetoelastic coupling is taken into account at the linear level only of
the strain tensor, by the coefficients, $B_{0}$ and $B_{1}$, that can be
derived either experimentally from magnetostriction constants Du Trémolet de
Lacheisserie (1993) or by ab initio calculations Fähnle  _et al._ (2002). It
should be stressed that the elastic constants, $C_{0}$ and $C_{1}$ can be
deduced from similar experiments, also.

Notwithstanding the presence of only the linear term of the strain, it does
describe non-trivial effects, that can be measured in real and numerical
experiments, as we have shown. In particular, the ultrafast switching of the
magnetization, induced by a damping-like STT in NiO, where both the duration
and intensity can be modulated by appropriate internal strains, thanks to the
magnetoelastic coupling.

It is, now, possible to set in perspective the limits of the present framework
and how they may be overcome.

The strain can be modulated, in return, by the induced magnetization process
by the same magnetoelastic coupling, so as to give rise to a modulated elastic
wave, whose properties, in extended systems can leave a characteristic
imprint, that would be interesting to look for.

It should be kept in mind, though, that, in real experiments, because
amplitude fluctuations of the electric pulses are inevitable, small damping
values are not favorable, since they may easily lead to overshoot, which is
amplified by the mechanical response. When Joule heating in the normal
underlying metal is taken into account, a shorter pulse with stronger current
intensity should be desirable, but a modulated undamped mechanical strain is
then produced. For the moment, a major conceptual issue is how to describe, in
a consistent way, the mechanical damping effects, that are induced by a
differential thermal stress, produced by such a Joule heating. We hope to
return to this issue in future work. Hints may be gleaned from our previous
work on the stochastic dynamics of the magnetization Tranchida  _et al._
(2018).

For the case of undamped dynamics we have developed the mathematical and
computational framework, that consistently takes into account elastic and
magnetic degrees of freedom and that provides a starting point for designing
future AFs devices for spintronics applications that combine both mechanical
and magnetic responses and their interactions. And we have, also, shown how to
take into account a particular class of damping terms, those that can describe
spin–transfer–torque sources.

It is this framework that is the starting point for taking into account the
stochastic effects mentioned above.

## Acknowledgments

All authors contributed equally to the manuscript. T.N. acknowledges financial
support through a joint doctoral fellowship “CEA–Région Centre” under grant
number 201600110872.

##  Appendix A A model for precession through anticommuting variables

We review some issues when dealing with classical spin variables and explain
some limitations which are met by using the framework of commuting variables
Brink  _et al._ (1976); *brinkLagrangianFormulationClassical1977. A common
starting point to build precessional models is choosing the appropriate
representation of the ${\mathfrak{su}}(2)$ algebra.

In the Heisenberg picture the Ehrenfest theorem implies that the equations of
motion for the average of the spin operator $\hat{\bm{S}}$ components read

$\imath\hbar\frac{d\langle\hat{\bm{S}}\rangle}{dt}=\langle[\hat{\bm{S}},{\cal
H}]\rangle$ (42)

where ${\cal H}$ is the Zeeman Hamiltonian operator for an external magnetic
field $\bm{H}$ such that

${\cal H}=-g\mu_{0}\mu_{B}\hat{\bm{S}}\cdot\bm{H}.$ (43)

From the commutation relations
$[\hat{S}_{i},\hat{S}_{j}]=\imath\epsilon_{ijk}\hat{S}_{k}$, one can quite
straightforwardly show that

$\frac{d\langle\hat{\bm{S}}\rangle}{dt}=\frac{g\mu_{B}\mu_{0}}{\hbar}\langle\hat{\bm{S}}\rangle\times\bm{H},$
(44)

which corresponds to the well-known Larmor precession of the expectation value
of the spin in an external magnetic field $\bm{H}$.

Interestingly, a different–but equivalent–approach can be followed, by
considering the Majorana representation, in terms of anticommuting variables,
$\xi_{I}$, $\xi_{I}\xi_{J}+\xi_{J}\xi_{I}=0$ Alicea (2012) of the vector
$\bm{S}$:

$S_{I}=-\frac{\imath}{2}\epsilon_{IJK}\xi_{J}\xi_{K},$ (45)

which can easily be shown to commute, as $S_{I}S_{J}=0$. This last relation
is, of course, not satisfactory. It can be easily avoided, however, by
imposing that the $\xi_{I}$ satisfy the anticommutation relations
$\\{\xi_{I},\xi_{J}\\}=\delta^{ij}$–that they generate a Clifford algebra (up
to a constant normalization) Lounesto (2003).

These relations imply that the $S_{I}$ define a spin–$\frac{1}{2}$
representation of the $\mathfrak{su}(2)$ algebra. Higher spin representations
can be defined by multilinear combinations Berezin and Marinov (1977), that
are relevant for describing magnetic properties of composite objects; since we
can work with the $S_{I}$ instead of the $\xi_{I}$, however, this complication
will not affect us here.

This representation highlights that the spin degrees of freedom, $S_{I}$, are,
in fact, “composite” objects and that the “fundamental degrees of freedom” are
the $\xi_{I}$. Therefore, it is useful to develop the description of the
dynamics, directly, in terms of the $\xi_{I}$ themselves. We shall recall the
salient features below.

The Poisson brackets for the anticommuting variables $\xi_{I}$ are related to
those of $S_{I}$, in order that the dynamics be, indeed, equivalent, in a way
that was set forth many years ago, through the construction of a corresponding
graded Poisson bracket, which generalizes Poisson brackets from manifolds to
super-manifolds Casalbuoni (1976a).

In terms of any functions of anti-commuting variables $\bm{\xi}$ this bracket
reads

$\left\\{f({\bm{\xi}}),g({\bm{\xi}})\right\\}_{\textrm{PB}}\equiv\frac{\imath}{\hbar}f({\bm{\xi}})\frac{\overleftarrow{\partial}}{\partial\xi_{K}}\frac{\overrightarrow{\partial}}{\partial\xi_{K}}g({\bm{\xi}}),$
(46)

with the corresponding definition of the left and right derivative of any
function of the anti-commuting variables. One can show that this bracket, also
known as the “antibracket” of any two functions on a flat supermanifold,
satisfies all necessary properties for a graded bracket, namely (graded)
Leibniz rule, (graded) anti-symmetry and (graded) Jacobi identity Gomis  _et
al._ (1995); *kanatchikovFieldTheoreticGeneralizations1997.

By taking this graded Poisson bracket for any two of these anti-commuting
variables, we get

$\\{\xi_{I},\xi_{J}\\}_{PB}=\frac{\imath}{\hbar}\delta_{IJ}.$ (47)

This implies that any two such variables are canonically conjugate, since
$\\{\xi_{I},-\imath\hbar\xi_{J}\\}_{PB}=\delta_{IJ}$ and
$\pi_{I}\equiv-\imath\hbar\xi_{I}$ defines the canonical conjugate Casalbuoni
(1976b); Berezin and Marinov (1977).

By using the Grassmann properties of ${\bm{\xi}}$, one proves that

$\\{S_{I},S_{J}\\}_{PB}=\frac{1}{\hbar}\epsilon_{IJK}S_{K}$ (48)

which is a consistency check that the $S_{I}$ as defined in (45) do realize a
representation of the rotation groupVarshalovich  _et al._ (1988).

In eq.(46), if the functions of the $\bm{\xi}$ are chosen to contain only
quadratic terms, then one can identify the previous bracket as a regular
Poisson bracket on a Riemannian manifold for the commuting variables $\bm{S}$

$\begin{array}[]{rcl}\displaystyle\\{f({\bm{S}}),g({\bm{S}})\\}_{PB}=\frac{\imath}{\hbar}\frac{\partial
f}{\partial
S_{I}}S_{I}\frac{\overleftarrow{\partial}}{\partial\xi_{K}}\frac{\overrightarrow{\partial}}{\partial\xi_{K}}S_{J}\frac{\partial
g}{\partial S_{J}}\\\ \displaystyle=\\{S_{I},S_{J}\\}_{PB}\frac{\partial
f}{\partial S_{I}}\frac{\partial g}{\partial
S_{J}}=\frac{1}{\hbar}\epsilon_{IJK}S_{K}\frac{\partial f}{\partial
S_{I}}\frac{\partial g}{\partial S_{J}},\end{array}$ (49)

which is, precisely, the “spinning part” of the bracket introduced by Yang and
Hirschfelder Yang and Hirschfelder (1980) for magnetized fluid dynamics, which
reads as

$\\{A,B\\}_{PB}\equiv\frac{\partial A}{\partial q_{I}}\frac{\partial
B}{\partial p_{I}}-\frac{\partial A}{\partial p_{I}}\frac{\partial B}{\partial
q_{I}}-\frac{1}{\hbar}\epsilon_{IJK}S_{I}\frac{\partial A}{\partial
S_{J}}\frac{\partial B}{\partial S_{K}}.$ (50)

Similar expressions have been found already by Casalbuoni a long time ago
Casalbuoni (1976b) without having attracted the attention they deserve.
Following the inverse path of canonical quantization, we use this graded
Poisson Bracket on any commuting quantity, which allows us to compute directly

$\\{\bm{S},{\cal H}\\}_{PB}=\frac{g\mu_{B}\mu_{0}}{\hbar}\bm{S}\times\bm{H}$
(51)

thereby highlighting that the description of $\bm{S}$ in terms of $\bm{\xi}$
is an equivalent description of the dynamics.

For any time-dependent commuting functions of ${\bm{S}}(t)$, we can use this
Poisson bracket to deduce a Liouville equation

$\frac{dF({\bm{S}}(t))}{dt}=\\{F({\bm{S}}(t)),{\cal H}\\}_{PB}.$ (52)

Consequently, the equations of motion for the variables ${\bm{\xi}}$ read much
more simple expressions as

$\frac{d{\bm{\xi}}(t)}{dt}=\\{{\bm{\xi}},{\cal
H}\\}_{PB}=\frac{\imath}{\hbar}\frac{\overrightarrow{\partial}{\cal
H}}{\partial{\bm{\xi}}},$ (53)

which are known to form a non-relativistic pseudoclassical mechanicsCasalbuoni
(1976b).

One reason why using the representation of $\bm{S}$ in terms of $\bm{\xi}$ is
useful is that it is intrinsically difficult to build a Lagrangian model for
the commuting spin variable ${\bm{S}}$, since its canonically conjugate
variable cannot be unambiguously identified. As is, by now, well known, the
conjugate of the dynamical variable ${\bm{\xi}}$ is proportional to itself
Casalbuoni (1976c). Therefore, the dynamics of spinning degrees of freedom can
be described, either through a vector of commuting variables on a curved
manifold, or by a vector of anti-commuting variables on a flat–though
non–Riemannian–manifold. Finally, it has been found that the Majorana-fermion
representation of $\frac{1}{2}$-spin operators is also a powerful tool to
straightforwardly compute spin-spin correlators Shnirman and Makhlin (2003);
*maoSpinDynamicsMajorana2003; *schadMajoranaRepresentationDissipative2015,
which represents an advantage for computing magnetic response functions in
many-body systems.

##  Appendix B Numerical accuracy of the splitting algorithms

The accuracy of the numerical schemes represented in table 1 depends on the
relative amplitude of the velocity terms
$(\dot{\varepsilon},\dot{\pi},\dot{S})$, and can be monitored by checking the
stability of equation (20) over time Hairer  _et al._ (2006); Blanes  _et al._
(2009).

Figure 6: (Color online) Single-site total energy as a function of both time
and variable timestep $Q$. Conditions of the simulation are expressed in
reduced units: $2C_{0}/\mu_{0}M_{s}^{2}=5.1\times 10^{5}$,
$2C_{1}/\mu_{0}M_{s}^{2}=3.5\times 10^{5}$,
$M_{1}V_{0}\mu_{0}M_{s}^{2}/2\hbar^{2}=1000$, ${\bm{\omega}}_{DC}=(0,0,2\pi)$
rad.GHz, $\pi_{11}(0)=1$, ${\bm{s}}(0)=(1,0,0)$. All the other parameters not
reported, included initial conditions are zero.

Figure 6 displays the total energy in one domain (given by the sum of eqs.(21)
and (22)) upon varying the numerical precision. The splitting algorithm
considered corresponds to the label $A=1$ in table 1. The numerical precision
is controlled by the “quality factor $Q$”, chosen at the beginning of the
simulation, which produces variable timesteps $\tau$ according to the relation
$\tau=Q/\|\bf{\omega}\|$. Simulation conditions produce a total energy equal
to $\pi_{11}^{2}(0)/4M_{1}$, which has to stay constant over time. We observe
that this is, indeed, the case, whatever the value of the quality factor $Q$.
We remark that, as $Q\to 0$, the variations about the average value of the
energy are suppressed, as should be expected.

Figure 7 displays some of the non-vanishing components of the strain tensor
over time (here $\varepsilon_{22}(t)=\varepsilon_{33}(t)$ and this third
component is not reported), when using different splitting algorithms, among
those displayed in Table 1. We observe less than $1\%$ of numerical relative
difference between the two algorithms on the strain and its conjugate
variables, and no difference (up to the machine precision) on the
magnetization with a “coarse” quality factor $Q$, which cannot be detected by
visual inspection on the figure. This difference falls to $0.1\%$ when the
quality factor is divided by $4$. The same procedure can be repeated for all
the splitting combinations in Table 1, the conclusions previously drawn apply,
also, for the magnetization, strain and its conjugate variable, depending on
the frequency of appearance of the splitted operator.

Figure 7: (Color online) Single-site strain components as a function of time
for various numerical schemes. Conditions of the simulation are identical to
those for figure 6 and the results are produced for $Q=0.0025$ only.

As expected, once a splitting algorithm is selected, the more frequently an
operator is evaluated, the smaller the error, though, of course, the time
required, also, grows. This opens the possibility to select optimally a proper
splitting algorithm, depending on the relative intensity of
$\dot{\varepsilon}_{ij},\dot{\pi}_{ij},{\dot{S}}_{i}$ over time.

## References

  * Radu  _et al._ (2011) I. Radu, K. Vahaplar, C. Stamm, T. Kachel, N. Pontius, H. Dürr, T. Ostler, J. Barker, R. Evans, R. Chantrell, _et al._ , Nature 472, 205 (2011).
  * Jungwirth  _et al._ (2016) T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich, Nature Nanotechnology 11, 231 (2016).
  * Železný  _et al._ (2018) J. Železný, P. Wadley, K. Olejník, A. Hoffmann, and H. Ohno, Nature Physics 14, 220 (2018).
  * Ahmad  _et al._ (2015) H. Ahmad, J. Atulasimha, and S. Bandyopadhyay, Scientific Reports 5, 18264 (2015).
  * Yan  _et al._ (2019) H. Yan, Z. Feng, S. Shang, X. Wang, Z. Hu, J. Wang, Z. Zhu, H. Wang, Z. Chen, H. Hua, W. Lu, J. Wang, P. Qin, H. Guo, X. Zhou, Z. Leng, Z. Liu, C. Jiang, M. Coey, and Z. Liu, Nature Nanotechnology 14, 131 (2019).
  * Fiebig (2005) M. Fiebig, Journal of Physics D: Applied Physics 38, R123 (2005).
  * Atulasimha and Bandyopadhyay (2010) J. Atulasimha and S. Bandyopadhyay, Applied Physics Letters 97, 173105 (2010).
  * Du Trémolet de Lacheisserie (1993) E. Du Trémolet de Lacheisserie,  _Magnetostriction: Theory and Applications of Magnetoelasticity_ (CRC Press, Boca Raton, 1993).
  * Maugin (1988) G. A. Maugin,  _Continuum Mechanics of Electromagnetic Solids_ , North Holland Series in Applied Mathematics and Mechanics No. 33 (North Holland, Amsterdam, 1988).
  * Kittel (1949) C. Kittel, Reviews of Modern Physics 21, 541 (1949).
  * Buiron  _et al._ (1999) N. Buiron, L. Hirsinger, and R. Billardon, Le Journal de Physique IV 09, Pr9 (1999).
  * Daniel and Galopin (2008) L. Daniel and N. Galopin, The European Physical Journal Applied Physics 42, 153 (2008).
  * Mei and Vernescu (2010) C. C. Mei and B. Vernescu,  _Homogenization Methods for Multiscale Mechanics_ (World Scientific, Singapore, 2010).
  * Casalbuoni (1976a) R. Casalbuoni, Il Nuovo Cimento A (1965-1970) 33, 115 (1976a).
  * Casalbuoni (1976b) R. Casalbuoni, Il Nuovo Cimento A (1965-1970) 33, 389 (1976b).
  * Berezin and Marinov (1977) F. Berezin and M. Marinov, Annals of Physics 104, 336 (1977).
  * Nelson and Chen (1994) D. F. Nelson and B. Chen, Physical Review B 50, 1023 (1994).
  * Cheng  _et al._ (2015) R. Cheng, M. W. Daniels, J.-G. Zhu, and D. Xiao, Physical Review B 91, 064423 (2015).
  * Nussle  _et al._ (2019) T. Nussle, P. Thibaudeau, and S. Nicolis, Journal of Magnetism and Magnetic Materials 469, 633 (2019), arXiv:1711.08062 .
  * Casimir (1945) H. B. G. Casimir, Reviews of Modern Physics 17, 343 (1945).
  * Landau  _et al._ (2008) L. D. Landau, E. M. Lifshitz, and L. D. Landau,  _Theory of Elasticity_ , 3rd ed., Course of Theoretical Physics No. Vol. 7 (Elsevier, Amsterdam, 2008).
  * Gairola (1978) B. K. D. Gairola, Physica Status Solidi (b) 85, 577 (1978).
  * Eringen (2004) A. C. Eringen, ed.,  _Nonlocal Continuum Field Theories_ (Springer New York, New York, NY, 2004).
  * Toupin and Bernstein (1961) R. A. Toupin and B. Bernstein, The Journal of the Acoustical Society of America 33, 216 (1961).
  * Callen and Callen (1963) E. R. Callen and H. B. Callen, Physical Review 129, 578 (1963).
  * Brown (1966) W. F. Brown,  _Magnetoelastic Interactions_, english ed., edited by C. Truesdell, R. Aris, L. Collatz, G. Fichera, P. Germain, J. Keller, M. M. Schiffer, and A. Seeger, Springer Tracts in Natural Philosophy, Vol. 9 (Springer Berlin Heidelberg, Berlin, Heidelberg, 1966).
  * Rosen (1972) G. Rosen, American Journal of Physics 40, 683 (1972).
  * Tiersten (1964) H. F. Tiersten, Journal of Mathematical Physics 5, 1298 (1964).
  * Akhiezer  _et al._ (1968) A. Akhiezer, V. Bar’yakhtar, and S. Peletminskii,  _Spin Waves_ , North Holland Series in Low Temperature Physics, Vol. 1 (Amsterdam, 1968).
  * Dever (1972) D. Dever, Journal of Applied Physics 43, 3293 (1972).
  * Wilczek (2009) F. Wilczek, Nature Physics 5, 614 (2009).
  * Alicea (2012) J. Alicea, Reports on Progress in Physics 75, 076501 (2012).
  * Lounesto (2003) P. Lounesto,  _Clifford Algebras and Spinors_ , 2nd ed., London Mathematical Society Lecture Note Series No. 286 (Cambridge Univ. Press, Cambridge, 2003).
  * Varshalovich  _et al._ (1988) D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii,  _Quantum Theory of Angular Momentum_ (World Scientific, 1988).
  * Evans  _et al._ (2014) R. F. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. Ellis, and R. W. Chantrell, Journal of Physics: Condensed Matter 26, 103202 (2014).
  * Eriksson  _et al._ (2017) O. Eriksson, A. Bergman, L. Bergqvist, and J. Hellsvik,  _Atomistic Spin Dynamics: Foundations and Applications_ , first edition ed. (Oxford University Press, Oxford, 2017).
  * van Vleck (1937) J. H. van Vleck, Physical Review 52, 1178 (1937).
  * Néel (1954) L. Néel, Journal de Physique et le Radium 15, 225 (1954).
  * Lee (1955) E. W. Lee, Reports on Progress in Physics 18, 184 (1955).
  * Sabiryanov and Jaswal (1999) R. F. Sabiryanov and S. S. Jaswal, Physical Review Letters 83, 2062 (1999).
  * Kastrup (1983) H. A. Kastrup, Physics Reports 101, 1 (1983).
  * Kanatchikov (1993) I. V. Kanatchikov, arXiv:9312162 (hep-th) (1993).
  * Arnowitt  _et al._ (1960) R. Arnowitt, S. Deser, and C. W. Misner, Physical Review 117, 1595 (1960).
  * Hojman  _et al._ (1976) S. A. Hojman, K. Kuchař, and C. Teitelboim, Annals of Physics 96, 88 (1976).
  * Hill (1970) R. Hill, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 314, 457 (1970).
  * Kijowski and Szczyrba (1976) J. Kijowski and W. Szczyrba, Communications in Mathematical Physics 46, 183 (1976).
  * Szczyrba (1976) W. Szczyrba, Communications in Mathematical Physics 51, 163 (1976).
  * Marsden  _et al._ (1986) J. E. Marsden, R. Montgomery, P. J. Morrison, and W. B. Thompson, Annals of Physics 169, 29 (1986).
  * Yang and Hirschfelder (1980) K.-H. Yang and J. O. Hirschfelder, Physical Review A 22, 1814 (1980).
  * Skomski (2008) R. Skomski,  _Simple Models of Magnetism_ , Oxford Graduate Texts (Oxford University Press, Oxford, New York, 2008).
  * Blanes  _et al._ (2009) S. Blanes, F. Casas, J. A. Oteo, and J. Ros, Physics Reports 470, 151 (2009).
  * Hairer  _et al._ (2006) E. Hairer, C. Lubich, and G. Wanner,  _Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations_ , 2nd ed., Springer Series in Computational Mathematics (Springer-Verlag, Berlin Heidelberg, 2006).
  * Batcho and Schlick (2001) P. F. Batcho and T. Schlick, The Journal of Chemical Physics 115, 4019 (2001).
  * Beaujouan  _et al._ (2012) D. Beaujouan, P. Thibaudeau, and C. Barreteau, Physical Review B 86, 174409 (2012).
  * Rodrigues (1840) O. Rodrigues, Journal de Mathématiques Pures et Appliquées 5, 380 (1840).
  * Honerkamp and Römer (1993) J. Honerkamp and H. Römer,  _Theoretical Physics: A Classical Approach_ (Springer, Berlin; New York, 1993).
  * Thibaudeau and Beaujouan (2012) P. Thibaudeau and D. Beaujouan, Physica A: Statistical Mechanics and its Applications 391, 1963 (2012).
  * Omelyan  _et al._ (2003) I. P. Omelyan, I. M. Mryglod, and R. Folk, Computer Physics Communications 151, 272 (2003).
  * Ogden (1997) R. W. Ogden,  _Non-Linear Elastic Deformations_ , dover ed ed., Dover Books on Physics (Dover, Mineola, NY, 1997).
  * Fähnle  _et al._ (2002) M. Fähnle, M. Komelj, R. Q. Wu, and G. Y. Guo, Physical Review B 65, 144436 (2002).
  * Tranchida  _et al._ (2018) J. Tranchida, P. Thibaudeau, and S. Nicolis, Physical Review E 98, 042101 (2018).
  * Brink  _et al._ (1976) L. Brink, P. Di Vecchia, and P. Howe, Physics Letters B 65, 471 (1976).
  * Brink  _et al._ (1977) L. Brink, P. Di Vecchia, and P. Howe, Nuclear Physics B 118, 76 (1977).
  * Gomis  _et al._ (1995) J. Gomis, J. París, and S. Samuel, Physics Reports 259, 1 (1995).
  * Kanatchikov (1997) I. V. Kanatchikov, Reports on Mathematical Physics 40, 225 (1997).
  * Casalbuoni (1976c) R. Casalbuoni, Physics Letters B 62, 49 (1976c).
  * Shnirman and Makhlin (2003) A. Shnirman and Y. Makhlin, Physical Review Letters 91, 207204 (2003).
  * Mao  _et al._ (2003) W. Mao, P. Coleman, C. Hooley, and D. Langreth, Physical Review Letters 91, 207203 (2003).
  * Schad  _et al._ (2015) P. Schad, Y. Makhlin, B. N. Narozhny, G. Schön, and A. Shnirman, Annals of Physics 361, 401 (2015).