Dynamic magnetostriction for antiferromagnets

TL;DR

This paper develops a Hamiltonian framework to analyze the dynamics of antiferromagnetic moments interacting with elastic media, revealing how magnetoelastic coupling influences switching behavior and magnetic oscillations.

Contribution

It introduces a symplectic geometric integration scheme for magnetoelastic dynamics in antiferromagnets, incorporating constraints from anticommuting variables without directly handling them.

Findings

Magnetoelastic coupling doubles spin accumulation.

Mechanical oscillations affect Nél order parameter dynamics.

Stress-induced strain slightly accelerates magnetization switching.

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\'eel 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.

Tables

Table 1. Table 1: Decomposition table of symplectic integrators

Equations

$$\frac{\mathcal{H}}{V_0}=\frac{{\mathcal{H}}_{\mathrm{mech}}}{V_0}+\frac{{\mathcal{H}}_{\mathrm{ext}}}{V_0}$$

$$\frac{{\mathcal{H}}_{\text{mech}}}{V_0}=\frac{1}{2}\sum _{i=1}^N{C}_{IJKL}^i{ϵ}_{IJ}^i(t){ϵ}_{KL}^i(t),$$

$$C_{IJKL}=C_0{\delta }_{IJ}{\delta }_{KL}+C_1\left({\delta }_{IK}{\delta }_{JL}+{\delta }_{IL}{\delta }_{JK}\right).$$

$$\frac{{\mathcal{H}}_{\text{ext}}}{V_0}=-\sum _{i=1}^N{\sigma }_{IJ}^{\text{ext}}{ϵ}_{JI}^i.$$

$${\mathcal{H}}_{\text{ME}}$$

$$+B_{IJKLMN}^{(2)i,j}{ϵ}_{IJ}^i(t){ϵ}_{KL}^j(t)S_M^i(t)S_N^j(t)+\dots$$

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

$$C_{IJLK}^{\text{eff},i,j}\equiv C_{IJLK}^{ij}+\frac{1}{V_0}{B}_{IJKLMN}^{(2)i,j}{S}_M^i(t)S_N^j(t).$$

$$S_I^i\equiv -\frac{ı}{2}{ϵ}_{IJK}{\xi }_J^i{\xi }_K^i.$$

$$S_I^{\text{eff}}=⟨S_I^i{⟩}_i$$

$${\mathcal{H}}_{\text{ME}}=B_{IJKL}^{(1)}{ϵ}_{IJ}(t)S_K(t)S_L(t).$$

$$B_{IJKL}^{(1)}=B_0^{(1)}{\delta }_{IJ}{\delta }_{KL}+B_1^{(1)}\left({\delta }_{IK}{\delta }_{JL}+{\delta }_{IL}{\delta }_{JK}\right).$$

$${\pi }_{IJ}\equiv \frac{\partial \mathcal{L}}{\partial {\dot{\epsilon }}_{IJ}},$$

$$\mathcal{H}(\mathbf{\epsilon },\mathbf{\pi })={\pi }_{IJ}{\dot{\epsilon }}_{IJ}-\mathcal{L},$$

$$\{A,B{\}}_{PB}=\frac{\partial A}{\partial {\epsilon }_{IJ}}\frac{\partial B}{\partial {\pi }_{IJ}}-\frac{\partial A}{\partial {\pi }_{IJ}}\frac{\partial B}{\partial {\epsilon }_{IJ}}.$$

$$\{{\epsilon }_{IJ},{\pi }_{KL}{\}}_{PB}$$

$$\{{\epsilon }_{IJ},{\epsilon }_{KL}{\}}_{PB}$$

$$\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.$$

$$\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}$$

$$\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}$$

$$\frac{\partial {\dot{\epsilon }}_{IJ}}{\partial {\epsilon }_{IJ}}+\frac{\partial {\dot{\pi }}_{IJ}}{\partial {\pi }_{IJ}}+\frac{\partial {\dot{S}}_K}{\partial S_K}=0.$$

$$U=$$

$$-\hslash {\omega }_I{S}_I,$$

$${\mathcal{H}}_{\mathrm{kinetic}}=\frac{1}{2}{\pi }_{IJ}{M}_{IJKL}^{-1}{\pi }_{KL}$$

$$M_{IJKL}=M_0{\delta }_{IJ}{\delta }_{KL}+M_1\left({\delta }_{IK}{\delta }_{JL}+{\delta }_{IL}{\delta }_{JK}\right).$$

$$\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}$$

$$\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.$$

$$\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}$$

$$\dot{f}$$

$$\equiv \mathcal{L}f.$$

$$\mathcal{L}={\dot{\epsilon }}_{IJ}\frac{\partial }{\partial {\epsilon }_{IJ}}+{\dot{\pi }}_{IJ}\frac{\partial }{\partial {\pi }_{IJ}}+{\dot{S}}_I\frac{\partial }{\partial S_I}.$$

$$f(\mathbf{\epsilon }(t),\mathbf{\pi }(t),\mathbf{S}(t))=e^{\mathcal{L}t}f(\mathbf{\epsilon }(0),\mathbf{\pi }(0),\mathbf{S}(0)).$$

Full text

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,

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

The total Hamiltonian can, thus, be written as

[TABLE]

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.

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

and the equations of motion become

[TABLE]

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

[TABLE]

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

[TABLE]

The last line defines the (total) Liouville operator

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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}$.

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.

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 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.

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.