# Hamiltonian and Lagrangian Dynamical Matrix Approaches Applied to Magnetic Nanostructures

Authors: Roberto Zivieri; Giancarlo Consolo
Journal: Advances in Condensed Matter Physics
 (2012)

---

## Abstract

Two micromagnetic tools to study the spin dynamics are reviewed. Both approaches are based upon the so-called dynamical matrix method, a hybrid micromagnetic framework used to investigate the spin-wave normal modes of confined magnetic systems. The approach which was formulated first is the Hamiltonian-based dynamical matrix method. This method, used to investigate dynamic magnetic properties of conservative systems, was originally developed for studying spin excitations in isolated magnetic nanoparticles and it has been recently generalized to study the dynamics of periodic magnetic nanoparticles. The other one, the Lagrangian-based dynamical matrix method, was formulated as an extension of the previous one in order to include also dissipative effects. Such dissipative phenomena are associated not only to intrinsic but also to extrinsic damping caused by injection of a spin current in the form of spin-transfer torque. This method is very accurate in identifying spin modes that become unstable under the action of a spin current. The analytical development of the system of the linearized equations of motion leads to a complex generalized Hermitian eigenvalue problem in the Hamiltonian dynamical matrix method and to a non-Hermitian one in the Lagrangian approach. In both cases, such systems have to be solved numerically.


---

## Body

## 1. Introduction

In these last years, great attention has been given to the study of magnetization dynamics in laterally confined magnetic systems. It is well know that spin excitations are quantized due to the lateral confinement. The oscillations are the so-called normal modes, which represent a pattern of motion given by all the parts of the system oscillating sinusoidally with the same frequency and with the same phase relation. In this last decade, analytical models have given important contributions to understand the frequency spectrum of normal modes for different ground-state magnetizations [1–11]. However, some limitations due to the assumptions made for the determination of the equilibrium state, the boundary conditions, and the calculation of the energy contributions to normal modes dynamics are still present.On the other hand, a lot of efforts have been devoted to develop micromagnetic codes having the aim of calculating very precisely different ground states of nanometric particles [12]. Due to their accuracy, the developed micromagnetic methods have contributed to give additional information about the spin dynamics. The first micromagnetic calculations were typically based upon codes developed to calculate in the first place the ground state of a given magnetic particle. Then the time evolution of the average magnetization of a particle could be obtained and, from a subsequent postprocessing of these data (mainly using the Fourier transform of the magnetization), information could be extracted about mode frequencies and spatial profiles [13]. In the simplest application of the method, the limit of these calculations was the observation of modes with nonzero magnetization only. More recently, a micromagnetic method was extended to the detailed calculation of eigenfrequencies and eigenvectors under the effect of an oscillatory in-plane small magnetic field [14]. Another recent micromagnetic method was also developed to study the quantized spin excitations in laterally confined systems [15]. This is the so-called Hamiltonian-based dynamical matrix method (HDMM).The HDMM was first formulated for isolated magnetic nanoelements and then it has been generalized to the case of interacting nanoelements [16]. It is the prototype of finite-difference methods and represents an eigenvalue/eigenvector problem. The scope of this method is to find the frequencies and profiles of the spin modes which are associated to the eigenvalues and the eigenvectors, respectively, of a dynamical matrix. It can be considered the analogous of the dynamical matrix formalism used to find atomic vibrations (phonons) in crystalline solids. The dynamical matrix contains the second derivatives of the density energy coming from a second order expansion of the density energy around the equilibrium. This method was already used to study the spin excitations in magnetic multilayers with ferro- or anti-ferromagnetic coupling [17, 18]. Due to the translational invariance, the number of independent dynamic variables was reduced to twice the number of the layers. The calculated second derivatives are evaluated at equilibrium. The eigenvalue/eigenvector problem can be set as a complex generalized Hermitian eigenvalue problem. The method presents several advantages: a single calculation yields the frequencies and eigenvectors of all modes of any symmetry, it is applicable to a particle of any shape (within the nanometric range), and the computation time is affordable. This means that by means of this micromagnetic approach, it is possible to determine, after a single iteration, the frequencies and the profiles of all spin-modes, independently of the ground-state magnetization (e.g., vortex state, vortex in the presence of an external magnetic field, onion state, quasi-saturated state). The main restriction of the method is its applicability to confined magnetic systems whose spin dynamics is assumed purely precessional with no dissipative effects. Of course, this is true only in a first approximation, since in real magnetic systems the intrinsic damping process plays an important rule. In order to select the representative modes of the spectrum and to compare them with the ones observed by means of the experimental techniques, the differential scattering cross-section has to be evaluated both for noninteracting and interacting magnetic particles.The first applications of the HDMM were on chains of dipolarly interacting rectangular dots representing a one-dimensional array [19] and of two-dimensional (2D) arrays formed by circular nanometric disks [20]. Very recently, this method was applied to study the collective mode dynamics in arrays of holes embedded into a thin ferromagnetic film. This calculation was done by including in the energy density computation also the exchange interaction between micromagnetic cells belonging to two adjacent primitive cells [21].In order to overcome the above-mentioned restrictions of the HDMM, Consolo et al. formulated very recently the so-called Lagrangian-based dynamical matrix method (LDMM) [22]. Such a method explicitly takes into account the intrinsic “positive” Gilbert damping and the current-induced spin-transfer-torque “negative” dissipation. Since the magnetic system so obtained is no more conservative, a Lagrangian formalism is necessary. Unlike the HDMM, the LDMM cannot be cast as a classical (nongeneralized) eigenvalue/eigenvector problem, but it has to be formulated as a complex generalized non-Hermitian eigenvalue problem. The first application of the method was done on a magnetic nanopillar stack of circular cross-section subject to an external magnetic field directed a few degree away from the normal to the plane [22]. The analysis was then extended to the case of external magnetic fields of variable intensity and orientation with respect to the plane of the nanopillar [23].It is important to notice that both formalisms have been developed up to now in the linear approximation, namely, considering small angular deviations of the magnetization from equilibrium so that each spin excitation is a normal mode of the system.The reference frame used in the micromagnetic calculations performed both by means of HDMM and LDMM is illustrated in Figure1. The z-axis is along the normal to the particle and the x-y plane lies on the particle plane. According to this reference frame, the configuration of the vector M→, representative of the magnetic dipole momentum, is identified through the polar angles, θ and ϕ, and the intrinsic rotation ψ. As it will appear clear in the section devoted to the Lagrangian approach, this latter angle does not enter in the equation of motion, being the corresponding Lagrange equation associated to a first integral of the motion (the conservation of the angular momentum). So that, as it is expected, because the modulus of the magnetization vector is preserved in time, the dynamics can be described through two degrees of freedom only (typically θ and ϕ).Figure 1
Reference frame used in the theory.If we assume that the magnetization vectorM→ is placed along the generic direction given by the unitary vector e^3′, it can be expressed in Cartesian coordinates by
(1)M→(t)=Ms(sin⁡θ(t)cos⁡ϕ(t),sin⁡θ(t)sin⁡ϕ(t),cos⁡θ(t)),
being Ms the saturation magnetization value (the modulus of the magnetization vector) and the time dependence of the angles θ and ϕ expressed as
(2)θ(t)=Θ+δθ(t),ϕ(t)=Φ+δϕ(t).
In (2), we indicate by (Θ, Φ) the static (equilibrium) orientation of the magnetization obtained by solving the stationary problem (M→.=0) for whatever effective field and magnetization distributions. Instead, δθ(t) and δϕ(t) are the small polar and azimuthal deviations from equilibrium, respectively.In the micromagnetic calculations, the magnetic system is subdivided into rectangular cells. Each micromagnetic cell is identified by a single indexk(j) that varies from 1 to N, where N denotes the number of cells. Hence, Mk is the magnetization in the kth cell and r→ij=r→i-r→j is the in-plane distance between the ith cell and the jth cell. The index has been assigned so that the first line of the rectangular matrix (X×Y) corresponds to k=1,…,Y and the second line to k=Y+1,…,2Y, and so on and so forth. X(Y) is the number of cells along x(y), while Z is the number of cells along z for a sample of thickness equal to d. We define for each cell the reduced magnetization m→k=M→k/Ms. Hence, in a polar reference frame for each cell
(3)m→k(t)=(sin⁡θk(t)cos⁡ϕk(t),sin⁡θk(t)sin⁡ϕk(t),cos⁡θk(t)),
where ϕk is the azimuthal angle and θk is the polar angle of the magnetization. The total energy density of the system, obtained by dividing the total energy by the volume of the cell, is a function of ϕk and θk : E=E(θk,ϕk) where k varies from 1 to N.In the following sections we describe in detail the two micromagnetic methods. We do not illustrate the applications of these methods to magnetic nanoparticles, because this aspect is not the purpose of this paper. We give a summary of the paper. In Section2, we outline the formalism at the basis of the HDMM. First, we introduce the different contributions to the energy density and then we derive the system of linear and homogeneous equations of motion for the isolated particle from the Hamilton equations. Finally, we present the generalization of the HDMM to the case of interacting elements. In Section 3, the formalism at the basis of the LDMM is presented. First, the Lagrangian equations for a macrospin system are derived. Then a generalization is given by considering N interacting momenta in an isolated magnetic element.


## 2. HDMM

This section deals with the review of the HDMM formalism. As stated above, the HDMM is a micromagnetic approach that can be applied only to fully conservative systems which are supposed to have, in a first approximation, a purely (undamped) precessional motion of the magnetization about the effective field. It is a finite-difference micromagnetic method developed in the linear regime of spin dynamics by considering small deviations from the equilibrium magnetization. In the first subsection, the different contributions of the energy density entering into the dynamical matrix are calculated and the equations of motion for an isolated magnetic element are cast in the form of a linear and homogeneous system that can be solved as a complex generalized Hermitian eigenvalue problem. In the last subsection, it is shown the generalization of the HDMM to interacting magnetic nanoparticles by including into the dynamical matrix the Bloch condition.


### 2.1. Energy Density

First, we give the explicit expressions of the different interactions entering into the total micromagnetic energy densityE of a given confined magnetic system simulated by using the HDMM: Zeeman energy, exchange energy, demagnetizing energy, and anisotropy energy, respectively [24]. The energy density is defined as E=E~/V where E~ is the energy of the system and V its volume. In the presence of an external magnetic field H→, the Zeeman energy density can be written in the form
(4)Eext=-μ0MsH→⋅∑k=1Nm→i,
being μ0 the vacuum permeability. In micromagnetic theory, the exchange energy can be expressed as a volume integral of the form
(5)E~exch=A∫part∑j=13(∇→mj)2dV,
where the subscript “part” denotes the volume of a general magnetic particle, A is the exchange stiffness constant, and ∇→ is the gradient applied to a given component of the magnetization. In this case the exchange contribution is independent of z. Using the first-neighbours model, the exchange energy density can be written as follows:
(6)Eexch=A∑k=1N∑n=141-m→k⋅m→nakn2,
where the variable akn is the distance between the centers of two adjacent cells of index k and n, respectively, k varies over all micromagnetic cells, and the sum over n ranges over the neighbours of the k-th cell. If the micromagnetic cells k are situated at the edges, one must impose boundary conditions. The cells on the edges interact with an external row of cells that have the same fixed magnetization, in this case the corresponding term in the sum must be weighted twice.In order to calculate the demagnetizing energy density, we have followed the method of the demagnetizing tensor. In the following, we use also the term dipolar in place of the term demagnetizing, because the higher-order terms of the expansion vanish in the practical cases examined. Generally, within the framework of the demagnetizing tensor method, the demagnetizing energy density can be written as follows:(7)Edmg=12μ0∑kjM→k⋅N__M→j=μ0Ms22∑kj(mxk,myk,mzk)×(NxxNxyNxzNyxNyyNyzNzxNzyNzz)(mxjmyjmzj).
This equation includes the self-energy and Nαβ=Nαβ(r→kj) with α,β=x,y,z are the elements of demagnetizing tensor. Each component of the demagnetizing tensor is related to the interaction between two rectangular surfaces S and S′. Under the assumption made for the calculation of the demagnetizing field for uniform magnetization, by using a version of Gauss’s theorem, the demagnetizing tensor can be written as
(8)N__(r→kj)=1V∫SkdS→∫Sj′dS→′|r→-r→′|,
where V=lc2d is the volume of the micromagnetic cell with lc the cell size and d is the cell height. Because of the four-fold C4 symmetry and since r→kj≡(x,y,z), all components can be expressed only as a function of Nxx(x,y,z) and Nxy(x,y,z) components with suitable permutations of the variables x,y,z.The magnetocrystalline uniaxial anisotropy can be labeled with the symbolEani. It is an energy density function, for a given micromagnetic cell, of the angle αk between the magnetization of the single cell m→k and the easy axis of generic direction given by the unit vector u^. We write
(9)Eani=∑k=1NK(1)sin2αk=∑k=1NK(1)(1-cos⁡2αk)=∑k=1NK(1)[1-(m→k⋅u^)2],
where K(1) is the first-order anisotropy uniaxial coefficient.As the dynamical matrix components are expressed in terms of the second derivatives of the energy density, it is necessary to calculate them from the above expressions. First, we calculate the second derivatives of the magnetization with respect to the polar and azimuthal angles of the given micromagnetic cell that represent the degrees of freedom of the system. Indeed, the second derivatives of the magnetization appear in the final expressions of the second derivatives of the energy density. In particular(10)∂2m→k∂δϕk2=(-sin⁡θkcos⁡ϕk,-sin⁡θksin⁡ϕk,0),∂2m→k∂δϕk∂δθk=(-cos⁡θksin⁡ϕk,cos⁡θkcos⁡ϕk,0),∂2m→k∂δθk2=(-sin⁡θkcos⁡ϕk,-sin⁡θksin⁡ϕk,-cos⁡θk).
The second derivatives of the energy density are calculated at equilibrium. For the sake of simplicity, in the following, the derivatives are calculated with respect to δθk and δθl implying that one or two of the two generic variables could be also δϕ.The second derivative of Zeeman energy density becomes(11)∂2Eext∂δθk∂δθl={-μ0MsH→⋅∂2m→k∂δθk∂δθll=k0l≠k.As outlined previously, for the calculation of the exchange contribution, the nearest-neighbour model is taken into account. It is useful to give also the expression of the first derivative due to some important manipulations that have to be performed. The first derivative with respect toδθk includes in the sum a term in which i=k and thus n≠k and also the other terms with n=k and with i one of the nearest neighbours. Thanks to a proper change of indices in the second term, the following equation is obtained:
(12)∂Eexch∂δθk=-A∑n=141akn2∂m→k∂δθk⋅m→n-A∑n=141akn2m→n⋅∂m→k∂δθk=-2A∑n=141akn2∂m→k∂δθk⋅m→n,
where the sum over n is made up over the nearest-neighbour micromagnetic cells of the k-th cell. In the special case of the adopted first neighbours model, the second derivatives are
(13)∂2Eexch∂δθk∂δθl={-2A∑n=141akn2∂2m→k∂δθk∂δθl⋅m→nl=k-2A1akl2∂m→k∂δθk⋅∂m→l∂δθll,k:nearest-neighbour0k≠landnotnearest-neighbour.
We now pass to the calculation of the derivatives of the demagnetizing energy density. Due to their rather complicated form, we give in the following the expression not only of the second derivatives of the demagnetizing energy density but also of the first derivatives with respect to the generic variable δθi. The first derivative of the demagnetizing energy density (7) is
(14)∂Edmg∂δθi=μ0Ms2[12∑k≠im→k⋅N__(k,i)∂m→i∂δθi+12∑j≠i∂m→i∂δθi∂Edmg∂δθiμ0Ms2⋅N__(k,i)m→j+12∂∂δθi(m→i⋅N^(i,i)m→i)∑k≠im→k⋅N__(k,i)∂m→i∂δθi+]=μ0Ms2∑k=1k≠iNm→k⋅N__(k,i)∂m→i∂δθi+m→i⋅N__(i,i)∂m→i∂δθi.
In the sum, the contribution of the terms with the same index k=j=i has been separated; moreover, we have taken into account that the tensor N__ fulfils N__(k,i)=N__(i,k) and is symmetric (Nαβ=Nβα).Thanks to the previous consideration, it is possible to write(15)∂m→k∂δθk⋅N__(k,i)m→i=m→i⋅N__(k,i)∂m→k∂δθk,
and therefore,
(16)∂∂δθi(m→i⋅N__(i,i)m→i)=∂m→i∂δθi⋅N__(i,i)m→i+m→i⋅N__(i,i)∂m→i∂δθi=2m→i⋅N__(i,i)∂m→i∂δθi.
The second derivative must take into account the two cases: i≠l,i=l(j=i), namely,
(17)∂2Edmg∂δθl∂δθi={Ms2(∑k=1Nm→k⋅N__(k,l)∂2m→l∂δθl∂δθiMs2+∂m→l∂δθl⋅N__(k,l)∂m→i∂δθi)l=iMs2∂m→l∂δθl⋅N__(k,l)∂m→i∂δθil≠i.
The k=i term resulting from the derivative of the second term (in the expression of the first derivative given above) has been included in the sum over k.The last step is the calculation of the term associated with the anisotropy energy density. The second derivative of the anisotropy energy density can be written as(18)∂2Eani∂δϕj∂δθk={-2K(1)[(∂m→k∂δϕj⋅u^)⋅(∂m→k∂δθk⋅u^)+(m→k⋅u^)-2K(1)(1)⋅(∂2m→k∂δϕj∂δθk⋅u^)]forδθ,δϕk=j-2K(1)[(∂m→k∂δϕk⋅u^)2+(m→k⋅u^)⋅(∂2m→k∂δϕk2⋅u^)]K(1)K(1)+(m→k⋅u^)forδϕ,δϕorδθ,δθk=j02K(1)K(1)+(m→k⋅u^)forδϕ,δϕorδθ,δθk≠j.


### 2.2. Equations of Motion for an Isolated Magnetic Particle

The derivatives calculated previously are included into the dynamic equations. Indeed, the equations of motion can be cast into a linear and homogeneous system in which the second derivatives of the energy density calculated at equilibrium appear explicitly. It is well know that the equation of motion for a magnetic spin system which undergoes a purely precessional motion is the Landau-Lifshitz equation [25], expressed as a torque equation involving the effective field and the magnetization itself. Since our aim is to find the energy density in a conservative system, we derive the equations of motion from the Hamilton equations.


#### 2.2.1. HDMM for a Macrospin System

The equation of motion for the magnetic systems under investigation will be derived by following semiclassical approach. Also, the model will be first derived by considering the so-called macrospin approximation, where the material is thought as uniformly magnetized and represented by a single dipole momentum.As known from classical mechanics, in the presence of fixed constraints and conservative sources, the system Hamiltonianℋ coincides with the total mechanical energy E~, namely, ℋ≡E~=T-U, where T is the kinetic energy and U is the potential expressed as the opposite of the potential energy V. By defining the Lagrangian variables of the problem with qn, where n=1,2,… is the number of degrees of freedom corresponding to the dynamic variables, and the corresponding conjugate momenta with pn, the Hamilton equations in the 2n canonical variables (qn, pn) take the form [15]
(19)∂qn∂t=∂H∂pn,∂pn∂t=-∂H∂qn.
The direction of the magnetic dipole moment of thekth cell is given by (3). For the specific case, the dynamic variables referred to the k-cell are the small deviations from equilibrium of the azimuthal and polar angles given by
(20)q1=δϕ,q2=δθ,
where 1 (2) labels the first (second) variable.To determine the conjugate momenta, the expression of the angular momentum is needed. We recall the relation between the angular momentuml→Ω and the magnetic momentum μ→ by referring to a a rigid body with a fixed point Ω (see Figure 1), namely,
(21)l→Ω=1γμ→=vcMsγm→,
where γ is the gyromagnetic ratio and vc the volume of the magnetic moment to the case. Indeed, since q1 represents a rotation about thez-axis of the magnetic dipole, its conjugate momentum p1 corresponds to the z component of the variation of the angular momentum, namely,(22a)p1=δlz=vcMsγδmz=-vcMsγsin⁡θδθ,
where m→=(sinθcosϕ,sinθsinϕ,cosθ) is the unit magnetization vector, normalized to the saturation magnetization Ms. Following an analogous argument, the other momentum p2 can be determined. Indeed, q2 is a rotation of the dipole moment about an axis of unitary vector ϕ^′=-sinϕe^1+cosϕe^2 in the x-yplane at an angle ϕ^′=ϕ+π/2 from the x-axis. Therefore, p2 corresponds to the projection of the variation of the angular momentum along the ϕ^′ vector, namely,(22b)p2=δl→Ω⋅ϕ′^=vcMsγ(δm→⋅ϕ^′)=-vcMsγ(∂m→∂ϕδϕ+∂m→∂θδθ)⋅ϕ^′=vcMsγsin⁡θδϕ.By substituting (20), (22a), and (22b) into the first Hamilton equation (cf. (19)), the following system of equations is obtained:
(23)δϕ.=-γvcMssin⁡θHδθ,δθ.=γvcMssin⁡θHδϕ,
where the dot notation stands for the time derivative and, at the right-hand side, the first derivatives of the energy with respect to the mechanical variables appear. By introducing the energy density E=E~/V (keeping in mind that ℋ=E~) and expanding it in a power Taylor expansion around the equilibrium up to the second order, it yields
(24)E=E0+12[Eϕϕ(δϕ)2+2Eϕθδϕδθ+Eθθ(δθ)2],
where E0 is the constant zero-order term that is inessential, the first-order terms vanish at equilibrium, and Eαβ represent the second derivatives calculated at equilibrium (Eαβ=∂2E/∂δα∂δβ with α,β=ϕ,θ). By using (23) and (24), we obtain
(25)δϕ.=-γMssin⁡θ[Eθϕδϕ+Eθθδθ],δθ.=γMssin⁡θ[Eϕϕδϕ+Eϕθδθ].
By inserting the time dependence in the form eiωt, where ω is the angular frequency of the given collective mode, the system of equations of motion reads
(26)-Eθϕsin⁡θδϕ-Eθθsin⁡θδθ-λ~δϕ=0,Eϕϕsin⁡θδϕ+Eϕθsin⁡θδθ-λ~δθ=0.
The linear and homogeneous system of equations expressed in (26) can be written as an eigenvalue problem
(27)C__v→=λ~v→,
with λ~=i(Ms/γ)ω the complex eigenvalues of the problem,
(28)C__=[-Eθϕsin⁡θ-Eθθsin⁡θEϕϕsin⁡θEϕθsin⁡θ],
being C__ a real, but not symmetric, matrix and v→=(δϕ,δθ)T.However, the system of motion equations (26) can be also recast as a complex generalized Hermitian eigenvalue problem
(29)A__v→=λB__v→,
where B__ is a Hessian matrix expressed by the second derivatives of the energy density at equilibrium. In particular
(30)B__=[EϕϕEϕθEθϕEθθ],A__=[0isin⁡θ-isin⁡θ0].
It should be noticed that the matrix A__ is Hermitian, whereas B__ is real and symmetric, so that all the corresponding eigenvalues λ=γ/Msω are real quantities.


#### 2.2.2. HDMM for an Isolated System Composed byN Interacting Magnetic Momenta

Equation (24) can be generalized to the case of N interacting magnetic momenta (where each momentum is identified with a micromagnetic cell) taking the form
(31)E=E0+12∑n=1N∑l=1N[Eϕnϕlδϕnδϕl+2Eϕnθlδϕnδθl+Eθnθlδθnδθl].
In (31), the total energy density is given by E=Eext+Eexch+Edmg+Eani where the different contributions are expressed in (4), (6), (7), and (9), respectively. By substituting (31) into (23), we get
(32)δϕk.=-γMssin⁡θk∑l=1N[Eθkϕlδϕl+Eθkθlδθl],δθk.=γMssin⁡θk∑l=1N[Eϕkϕlδϕl+Eϕkθlδθl].
By introducing the time dependence in the form eiωt, the system of equations of motion is composed by the following 2Nlinear and homogeneous equations for k=1⋯N:
(33)∑l=1N(-Eθkϕlsin⁡θk)δϕl+∑l=1N(-Eθkθlsin⁡θk)δθl-λ~δϕk=0,∑l=1N(Eϕkϕlsin⁡θk)δϕl+∑l=1N(Eϕkθlsin⁡θk)δθl-λ~δθk=0.
The unknown factors δϕl,δθl represent the eigenvectors of the problem and are expressed by the small angular deviation from the equilibrium position of the azimuthal (ϕl) and polar (θl) angles in the lth micromagnetic cell. The system above has a solution only if the determinant is zero. By suitable exchanges of rows (columns), the linear and homogeneous system of equations expressed in (33) can be written as an eigenvalue problem in analogy with the case of a macrospin system
(34)C__v→=λ~v→.
In (34), v→ is the set of the unknown factors representing the eigenvectors of the problem that take the form
(35)v→=(δϕ1,δθ1,δϕ2,δθ2,…,δϕN,δθN)T.C__ is the matrix whose elements are expressed as
(36)C2k-1,2l-1=-Eθkϕlsin⁡θkC2k-1,2l=-Eθkθlsin⁡θkC2k,2l-1=Eϕkϕlsin⁡θkC2k,2l=Eϕkθlsin⁡θk}k=1⋯N,l=1⋯N.
Analogously to the case of a macrospin system, the matrix C__ is real but not symmetric, also in this case the eigenvalues (and not only the eigenvectors) are complex. This matrix can be indeed seen as composed by two submatrices 2×2 for each pair of values (k,l). In the diagonal submatrices (k=l) the following relation is verified:
(37)C2k-1,2k-1=-C2k,2k.
For the elements of two different submatrices (k,l) and (l,k), that are not diagonal (k≠l), the following symmetries hold:
(38)sin⁡θkC2k-1,2l-1=-sin⁡θlC2l,2k,sin⁡θkC2k-1,2l=sin⁡θlC2l-1,2k,sin⁡θkC2k,2l-1=sin⁡θlC2l,2k-1,sin⁡θkC2k,2l=-sin⁡θlC2l-1,2k-1.
As for the macrospin approximation, the equation of motion can be recast as a complex generalized Hermitian eigenvalue problem
(39)A__v→=λB__v→,
where B__ is a Hessian matrix expressed by the second derivatives of the energy density at equilibrium. B__ is given by
(40)B2k-1,2l-1=EϕkϕlB2k-1,2l=EϕkθlB2k,2l-1=EθkϕlB2k,2l=Eθkθl}k=1⋯N,l=1⋯N,
where the matrix B__ is, again, real and symmetric. Moreover, since the static magnetization corresponds to a minimum of the energy and the matrix B__ is its Hessian, the matrix B__ is also positive defined. Instead, the matrix A__ has the following form:
(41)A__=[0isin⁡θ100⋯-isin⁡θ1000⋯000isin⁡θ2⋯00-isin⁡θ20⋯⋯⋯].
The matrix A__ is Hermitian. This allows us to solve the system as a complex generalized Hermitian eigenvalue problem which admits only real eigenvalues. To further reduce the computational time, it is possible to evaluate only some eigenvalues and eigenvectors that are in a specific range.Once the eigenvectorsv→ are obtained, the dynamic magnetization δm→k in the kth micromagnetic cell expressed in Cartesian coordinates and in unit of MS is given by
(42)δm→k=(-sin⁡θksin⁡ϕkδϕk+cos⁡θkcos⁡ϕkδθk,sisin⁡θkcos⁡ϕkδϕk+cos⁡θksin⁡ϕkδθk,-sin⁡θkδθk).
For each solution of the eigenvalue problem, the collection of all δm→k defines the mode profile. It must be remarked that δm→k is a complex vector, because δθk,δϕk are, in general, complex.


### 2.3. Equations of Motion for Interacting Magnetic Particles

Let us suppose to have a 2D periodic array of interacting nanodots characterized by the primitive vectorsa→1 and a→2; for example, for the specific case of a rectangular lattice their values are
(43)a→1=λxx^,a→2=λyy^,
where λx and λy represent the periodicity along x-axis and y-axis, respectively. The primitive vectors of the reciprocal lattice are b→1 and b→2(44)b→1=2π(a→2×a→1)×a→2(a→1×a→2)2,b→2=2π(a→1×a→2)×a→2(a→1×a→2)2.
For the special case of a rectangular lattice the vectors become
(45)b→1=2πλxx^,b→2=2πλyy^.
Due to the analogy with the Bloch wave (analogy, not equality, because the wave function has not a physical meaning, contrarily to the magnetization), it is possible to write the following periodicity rule valid for the dynamic magnetization:
(46)δm→(r→+R→)=eiK→⋅R→δm→(r→),
where R→ is a vector of the particle lattice given by
(47)R→=i1a→1+i2a→2i1,i2∈Z,i1,i2=-Ni2⋯Ni2-1,
and r→ can be confined into the first primitive cell centered in the first dot. The Bloch vector takes the following values:
(48)K→=n1N1b→1+n2N2b→2,ni∈Z,ni=-Ni2⋯Ni2-1.N1,N2∈ℕ indicate the number of primitive cells n in direction a→1 and a→2, respectively. In this scheme, both N1 and N2 are taken as even numbers. In order to confirm the hypothesis on the dynamic magnetization N1,N2 must be very large.If the magnetizations of different primitive cells and the different micromagnetic cells were independent, that is no periodicity rule were present, then one would have a dynamic system with variablesθkR→ and ϕkR→, where the k index changes inside the magnetic particle and R→ can assume the values indicated in (47). In this case, the linear and homogeneous system of 2N equations of motion is
(49)∑l,R→′(-EθkR→ϕlR→′sin⁡θkR→)δϕlR→′+∑l,R′→(-EθkR→θlR→′sin⁡θjkR→)δθlR→′-λδϕkR→=0,∑l,R→′(EϕkR→ϕlR→′sin⁡θkR→)δϕlR→′+∑l,R′→(EϕkR→θlR→′sin⁡θkR→)δθlR→′-λδθkR→=0,
where k=1⋯N and the sums over l and R→′ are on the same values. Instead, thanks to the Bloch condition expressed in (46), one can consider the equations only at R→=0; moreover, taking into account the same condition, the variables appearing for R→′≠0 can be replaced by using the same condition. Now, when rewriting the system, the index R→ is omitted when it has value equal to 0 or it is irrelevant. Owing to these considerations, the system given in (49) can be rewritten in the form [16]
(50)∑l=1N(-∑R→′EθkϕlR→′eiK→⋅R→′sin⁡θk)δϕl+∑l=1N(-∑R→′EθkθlR→′eiK→⋅R→′sin⁡θk)δθlδθl-λ~δϕk=0,∑l=1N(∑R→′EϕkϕlR→′eiK→⋅R→′sin⁡θk)δϕl+∑l=1N(∑R→′EϕkθlR→′eiK→⋅R→′sin⁡θk)δθlδθl-λ~δθk=0.
Equation (50) is similar to (33) by making the following replacement:
(51)Eαkβl⟶∑R→′eiK→⋅R→′Eαk0βlR→′,
and recalling that now the energy is referred to the whole system of particles. Like for the case of the isolated magnetic particle, also for the case of interacting magnetic particles the system of linear and homogeneous equations given in (50) can be written as an eigenvalue problem which in turn can be cast as a complex generalized Hermitian eigenvalue problem.The symmetry for the matrix elements that was valid for a single primitive cell now is not respected except forK→=0 or K→=G→/2 with G→ a translational reciprocal vector:
(52)∑R→′eiK→⋅R→′Eαk0→βlR→′=∑R→′eiK→⋅R→′EβlR→′αk0→=∑R→′eiK→⋅R→′Eβl0→αk-R→′=∑R→′e-iK→⋅R→″Eβl0→αkR→″≠∑R→′eiK→⋅R→′Eβl0→αkR→′,K→≠0→,K→≠G→2,
with R→″=-R→′. The primitive cell has at the centre a single dot that occupies only a part of it. The interdot exchange coupling is zero. Thanks to the last consideration and to the fact that derivatives of Zeeman, exchange, anisotropy, energy density are referred only to the cell of the first variable (αk0→) with αkR→=θkR→,ϕkR→ or at most to the nearest neighbour, all terms of the sum in (51) with R→′≠0 are zero. Hence, for these energy density terms, the equations are the same as those of the single particle case and the same occurs for their corresponding derivatives appearing in the equations of motion. The only energy density term that differs from the one obtained for the isolated element is the demagnetizing energy density. For a system of interacting nanoparticles, the demagnetizing energy density can be written as
(53)Edmg=12μ0∑R→,R′→,k,k′m→k(R→)N__(R→,R→′,k,k′)m→k′(R→′).
Due to its rather complex expression, it is useful to give the derivation also of the first derivative of Edmg like for the case of the isolated nanoelement.In order to calculate the first derivative, the properties of the demagnetizing tensor must be considered(54a)〈v|N|ω〉=〈ω|N+|v〉*=〈ω|N+|v〉=〈ω|Nt|v〉=〈ω|N|v〉,
since N is real and symmetric and(54b)N__(R→,R→′,k,k′)=N__(R→′-R→+r→k′-r→k)=N__(-R→′+R→-r→k′+r→k)=N__(R→′,R→,k′,k),
thanks to the inversion symmetry. Hence, the first derivative of the energy is(55)∂Edmg∂δαk(0→)μ0=μ0Ms22[∑R→′k′(R→′,k′)≠(0→,k)m→k′(R→′)⋅N__(R→′,0→,k′,k)∂m→k(0→)∂δαk(0→)μ0=μ0Ms22μ0+∑R→′k′(R→′,k′)≠(0→,k)∂m→k(0→)∂δαk(0→)⋅N__(0→,R→′,k′,k)m→k′(R→′)μ0=μ0Ms22μ0+∂∂δαk(0→)(m→k(0→)⋅N__(0→,0→,k,k)m→k(0→))∑R→′k′(R→′,k′)≠(0→,k)m→k′(R→′)⋅N__(R→′,0→,k′,k)∂m→k(0→)∂δαk(0→)]μ0=μ0Ms2[∑R→′k′(R→′,k′)≠(0→,k)m→k′(R→′)⋅N__(R→′,0→,k′,k)∂m→k(0→)∂δαk(0→)μ0=μ0Ms22μ0+m→k(0→)⋅N__(0→,0→,k,k)∂m→k(0→)∂δαk(0→)∑R→′k′(R→′,k′)≠(0→,k)m→k′(R→′)⋅N__(R→′,0→,k′,k)∂m→k(0→)∂δαk(0→)].
The second derivative takes the form
(56)∂2Edmg∂δαk(0→)∂δβk(R→)μ0={μ0MS2[∑R→′k′m→k′(R→′)⋅N__(R→′,0→,k′,k)μ0MS2μ0×∂2m→k(0→)∂δαk(0→)∂δβk(0→)μ0MS2μ0+∂m→k(0→)∂δβk(0→)μ0MS2μ0⋅N__(0→,0→,k,k)∂m→k(0→)∂δαk(0→)]R→=0→,l=kμ0Ms2∂m→l(0→)∂δβl(0→)⋅N__(R→,0→,l,k)∂m→k(0→)∂δαk(0→)(R→,l)≠(0→,k).
In the sum it is included the case in which (R→′,k′)=(0→,k) that is generated from the derivative of the second term in the expression of the first derivative given above. Now it is possible to calculate the terms that enter into the system of (50), starting with the one corresponding to l=k(57a)∑R→eiK→⋅R→Eαk0→βlR→=∑R→eiK→⋅R→Eαk0→βkR→∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2[∑R→′k′m→k′(0→)⋅N__(R→′,0→,k′,k)∑R→eiK→⋅R→Eαk0→βlR→×∂2m→k(0→)∂δαk(0→)∂δβk(0→)+∂m→k(0→)∂δβk(0→)∑R→eiK→⋅R→Eαk0→βlR→⋅N__(0→,0→,k,k)∂m→k(0→)∂δαk(0→)]∑R→eiK→⋅R→Eαk0→βlR→+μ0Ms2∑R→≠0→eiK→⋅R→∂m→k(0→)∂δβk(0→)∑R→eiK→⋅R→Eαk0→βlR→⋅N__(R→,0→,k,k)∂m→k(0→)∂δαk(0→)∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2∑R→[∑k′m→k′(0→)⋅N__(R→,0→,k′,k)∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2μ0Ms2×∂2m→k(0→)∂δαk(0→)∂δβk(0→)∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2μ0Ms2+eiK→⋅R→∂m→k(0→)∂δβk(0→)∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2μ0Ms2⋅N__(R→,0→,k,k)∂m→k(0→)∂δαk(0→)].
Equation (57a) is obtained by taking into account (46), by including the term corresponding to R→=0 into the sum performed over R→≠0 and by considering that both the static magnetization and the dynamic magnetization do not depend on R→. Indeed, second derivatives are calculated at equilibrium and the exponential eiK→·R→ appears on both the numerator and the denominator of the derivative.Whenl≠k, the term turns out to be
(57b)∑R→eiK→⋅R→Eαk0→βlR→μ0=μ0Ms2[∑R→eiK→⋅R→∂m→l(0→)∂δβl(0→)⋅N__(R→,0→,l,k)∂m→k(0→)∂δαk(0→)].Due to the properties of the demagnetizing tensor, the symmetry(58)∑R→′eiK→⋅R→′Eαk0→βlR→′=∑R→′eiK→⋅R→′Eαl0→βkR→′,
is fulfilled when l=k, but it is not fulfilled when l≠k.The formalism previously developed for interacting particles can be extended to a system of 2D antidots (ADs). In this case, it is necessary to add the exchange interaction between primitive cells [21]. In extended magnetic system like AD arrays, in addition to the usual nearest-neighbours exchange interaction between micromagnetic cells, the exchange contribution across the nearest-neighbours micromagnetic cells belonging to adjacent surface primitive cells must be taken into account. Hence, we recall the exchange energy density of (6)
(59)Eexch=A∑k∑n(1-m→k⋅m→n)akn2.
Here, the first sum runs over all the micromagnetic cells of the primitive cell and the second sum runs over the nearest neighbours of the kth micromagnetic cell. When the kth micromagnetic cell is on one of the edges (vertices) of the given primitive cell, the interaction with one (two) micromagnetic cell(s) belonging to the correct nearest primitive cell must be added.


## 2.1. Energy Density

First, we give the explicit expressions of the different interactions entering into the total micromagnetic energy densityE of a given confined magnetic system simulated by using the HDMM: Zeeman energy, exchange energy, demagnetizing energy, and anisotropy energy, respectively [24]. The energy density is defined as E=E~/V where E~ is the energy of the system and V its volume. In the presence of an external magnetic field H→, the Zeeman energy density can be written in the form
(4)Eext=-μ0MsH→⋅∑k=1Nm→i,
being μ0 the vacuum permeability. In micromagnetic theory, the exchange energy can be expressed as a volume integral of the form
(5)E~exch=A∫part∑j=13(∇→mj)2dV,
where the subscript “part” denotes the volume of a general magnetic particle, A is the exchange stiffness constant, and ∇→ is the gradient applied to a given component of the magnetization. In this case the exchange contribution is independent of z. Using the first-neighbours model, the exchange energy density can be written as follows:
(6)Eexch=A∑k=1N∑n=141-m→k⋅m→nakn2,
where the variable akn is the distance between the centers of two adjacent cells of index k and n, respectively, k varies over all micromagnetic cells, and the sum over n ranges over the neighbours of the k-th cell. If the micromagnetic cells k are situated at the edges, one must impose boundary conditions. The cells on the edges interact with an external row of cells that have the same fixed magnetization, in this case the corresponding term in the sum must be weighted twice.In order to calculate the demagnetizing energy density, we have followed the method of the demagnetizing tensor. In the following, we use also the term dipolar in place of the term demagnetizing, because the higher-order terms of the expansion vanish in the practical cases examined. Generally, within the framework of the demagnetizing tensor method, the demagnetizing energy density can be written as follows:(7)Edmg=12μ0∑kjM→k⋅N__M→j=μ0Ms22∑kj(mxk,myk,mzk)×(NxxNxyNxzNyxNyyNyzNzxNzyNzz)(mxjmyjmzj).
This equation includes the self-energy and Nαβ=Nαβ(r→kj) with α,β=x,y,z are the elements of demagnetizing tensor. Each component of the demagnetizing tensor is related to the interaction between two rectangular surfaces S and S′. Under the assumption made for the calculation of the demagnetizing field for uniform magnetization, by using a version of Gauss’s theorem, the demagnetizing tensor can be written as
(8)N__(r→kj)=1V∫SkdS→∫Sj′dS→′|r→-r→′|,
where V=lc2d is the volume of the micromagnetic cell with lc the cell size and d is the cell height. Because of the four-fold C4 symmetry and since r→kj≡(x,y,z), all components can be expressed only as a function of Nxx(x,y,z) and Nxy(x,y,z) components with suitable permutations of the variables x,y,z.The magnetocrystalline uniaxial anisotropy can be labeled with the symbolEani. It is an energy density function, for a given micromagnetic cell, of the angle αk between the magnetization of the single cell m→k and the easy axis of generic direction given by the unit vector u^. We write
(9)Eani=∑k=1NK(1)sin2αk=∑k=1NK(1)(1-cos⁡2αk)=∑k=1NK(1)[1-(m→k⋅u^)2],
where K(1) is the first-order anisotropy uniaxial coefficient.As the dynamical matrix components are expressed in terms of the second derivatives of the energy density, it is necessary to calculate them from the above expressions. First, we calculate the second derivatives of the magnetization with respect to the polar and azimuthal angles of the given micromagnetic cell that represent the degrees of freedom of the system. Indeed, the second derivatives of the magnetization appear in the final expressions of the second derivatives of the energy density. In particular(10)∂2m→k∂δϕk2=(-sin⁡θkcos⁡ϕk,-sin⁡θksin⁡ϕk,0),∂2m→k∂δϕk∂δθk=(-cos⁡θksin⁡ϕk,cos⁡θkcos⁡ϕk,0),∂2m→k∂δθk2=(-sin⁡θkcos⁡ϕk,-sin⁡θksin⁡ϕk,-cos⁡θk).
The second derivatives of the energy density are calculated at equilibrium. For the sake of simplicity, in the following, the derivatives are calculated with respect to δθk and δθl implying that one or two of the two generic variables could be also δϕ.The second derivative of Zeeman energy density becomes(11)∂2Eext∂δθk∂δθl={-μ0MsH→⋅∂2m→k∂δθk∂δθll=k0l≠k.As outlined previously, for the calculation of the exchange contribution, the nearest-neighbour model is taken into account. It is useful to give also the expression of the first derivative due to some important manipulations that have to be performed. The first derivative with respect toδθk includes in the sum a term in which i=k and thus n≠k and also the other terms with n=k and with i one of the nearest neighbours. Thanks to a proper change of indices in the second term, the following equation is obtained:
(12)∂Eexch∂δθk=-A∑n=141akn2∂m→k∂δθk⋅m→n-A∑n=141akn2m→n⋅∂m→k∂δθk=-2A∑n=141akn2∂m→k∂δθk⋅m→n,
where the sum over n is made up over the nearest-neighbour micromagnetic cells of the k-th cell. In the special case of the adopted first neighbours model, the second derivatives are
(13)∂2Eexch∂δθk∂δθl={-2A∑n=141akn2∂2m→k∂δθk∂δθl⋅m→nl=k-2A1akl2∂m→k∂δθk⋅∂m→l∂δθll,k:nearest-neighbour0k≠landnotnearest-neighbour.
We now pass to the calculation of the derivatives of the demagnetizing energy density. Due to their rather complicated form, we give in the following the expression not only of the second derivatives of the demagnetizing energy density but also of the first derivatives with respect to the generic variable δθi. The first derivative of the demagnetizing energy density (7) is
(14)∂Edmg∂δθi=μ0Ms2[12∑k≠im→k⋅N__(k,i)∂m→i∂δθi+12∑j≠i∂m→i∂δθi∂Edmg∂δθiμ0Ms2⋅N__(k,i)m→j+12∂∂δθi(m→i⋅N^(i,i)m→i)∑k≠im→k⋅N__(k,i)∂m→i∂δθi+]=μ0Ms2∑k=1k≠iNm→k⋅N__(k,i)∂m→i∂δθi+m→i⋅N__(i,i)∂m→i∂δθi.
In the sum, the contribution of the terms with the same index k=j=i has been separated; moreover, we have taken into account that the tensor N__ fulfils N__(k,i)=N__(i,k) and is symmetric (Nαβ=Nβα).Thanks to the previous consideration, it is possible to write(15)∂m→k∂δθk⋅N__(k,i)m→i=m→i⋅N__(k,i)∂m→k∂δθk,
and therefore,
(16)∂∂δθi(m→i⋅N__(i,i)m→i)=∂m→i∂δθi⋅N__(i,i)m→i+m→i⋅N__(i,i)∂m→i∂δθi=2m→i⋅N__(i,i)∂m→i∂δθi.
The second derivative must take into account the two cases: i≠l,i=l(j=i), namely,
(17)∂2Edmg∂δθl∂δθi={Ms2(∑k=1Nm→k⋅N__(k,l)∂2m→l∂δθl∂δθiMs2+∂m→l∂δθl⋅N__(k,l)∂m→i∂δθi)l=iMs2∂m→l∂δθl⋅N__(k,l)∂m→i∂δθil≠i.
The k=i term resulting from the derivative of the second term (in the expression of the first derivative given above) has been included in the sum over k.The last step is the calculation of the term associated with the anisotropy energy density. The second derivative of the anisotropy energy density can be written as(18)∂2Eani∂δϕj∂δθk={-2K(1)[(∂m→k∂δϕj⋅u^)⋅(∂m→k∂δθk⋅u^)+(m→k⋅u^)-2K(1)(1)⋅(∂2m→k∂δϕj∂δθk⋅u^)]forδθ,δϕk=j-2K(1)[(∂m→k∂δϕk⋅u^)2+(m→k⋅u^)⋅(∂2m→k∂δϕk2⋅u^)]K(1)K(1)+(m→k⋅u^)forδϕ,δϕorδθ,δθk=j02K(1)K(1)+(m→k⋅u^)forδϕ,δϕorδθ,δθk≠j.


## 2.2. Equations of Motion for an Isolated Magnetic Particle

The derivatives calculated previously are included into the dynamic equations. Indeed, the equations of motion can be cast into a linear and homogeneous system in which the second derivatives of the energy density calculated at equilibrium appear explicitly. It is well know that the equation of motion for a magnetic spin system which undergoes a purely precessional motion is the Landau-Lifshitz equation [25], expressed as a torque equation involving the effective field and the magnetization itself. Since our aim is to find the energy density in a conservative system, we derive the equations of motion from the Hamilton equations.


### 2.2.1. HDMM for a Macrospin System

The equation of motion for the magnetic systems under investigation will be derived by following semiclassical approach. Also, the model will be first derived by considering the so-called macrospin approximation, where the material is thought as uniformly magnetized and represented by a single dipole momentum.As known from classical mechanics, in the presence of fixed constraints and conservative sources, the system Hamiltonianℋ coincides with the total mechanical energy E~, namely, ℋ≡E~=T-U, where T is the kinetic energy and U is the potential expressed as the opposite of the potential energy V. By defining the Lagrangian variables of the problem with qn, where n=1,2,… is the number of degrees of freedom corresponding to the dynamic variables, and the corresponding conjugate momenta with pn, the Hamilton equations in the 2n canonical variables (qn, pn) take the form [15]
(19)∂qn∂t=∂H∂pn,∂pn∂t=-∂H∂qn.
The direction of the magnetic dipole moment of thekth cell is given by (3). For the specific case, the dynamic variables referred to the k-cell are the small deviations from equilibrium of the azimuthal and polar angles given by
(20)q1=δϕ,q2=δθ,
where 1 (2) labels the first (second) variable.To determine the conjugate momenta, the expression of the angular momentum is needed. We recall the relation between the angular momentuml→Ω and the magnetic momentum μ→ by referring to a a rigid body with a fixed point Ω (see Figure 1), namely,
(21)l→Ω=1γμ→=vcMsγm→,
where γ is the gyromagnetic ratio and vc the volume of the magnetic moment to the case. Indeed, since q1 represents a rotation about thez-axis of the magnetic dipole, its conjugate momentum p1 corresponds to the z component of the variation of the angular momentum, namely,(22a)p1=δlz=vcMsγδmz=-vcMsγsin⁡θδθ,
where m→=(sinθcosϕ,sinθsinϕ,cosθ) is the unit magnetization vector, normalized to the saturation magnetization Ms. Following an analogous argument, the other momentum p2 can be determined. Indeed, q2 is a rotation of the dipole moment about an axis of unitary vector ϕ^′=-sinϕe^1+cosϕe^2 in the x-yplane at an angle ϕ^′=ϕ+π/2 from the x-axis. Therefore, p2 corresponds to the projection of the variation of the angular momentum along the ϕ^′ vector, namely,(22b)p2=δl→Ω⋅ϕ′^=vcMsγ(δm→⋅ϕ^′)=-vcMsγ(∂m→∂ϕδϕ+∂m→∂θδθ)⋅ϕ^′=vcMsγsin⁡θδϕ.By substituting (20), (22a), and (22b) into the first Hamilton equation (cf. (19)), the following system of equations is obtained:
(23)δϕ.=-γvcMssin⁡θHδθ,δθ.=γvcMssin⁡θHδϕ,
where the dot notation stands for the time derivative and, at the right-hand side, the first derivatives of the energy with respect to the mechanical variables appear. By introducing the energy density E=E~/V (keeping in mind that ℋ=E~) and expanding it in a power Taylor expansion around the equilibrium up to the second order, it yields
(24)E=E0+12[Eϕϕ(δϕ)2+2Eϕθδϕδθ+Eθθ(δθ)2],
where E0 is the constant zero-order term that is inessential, the first-order terms vanish at equilibrium, and Eαβ represent the second derivatives calculated at equilibrium (Eαβ=∂2E/∂δα∂δβ with α,β=ϕ,θ). By using (23) and (24), we obtain
(25)δϕ.=-γMssin⁡θ[Eθϕδϕ+Eθθδθ],δθ.=γMssin⁡θ[Eϕϕδϕ+Eϕθδθ].
By inserting the time dependence in the form eiωt, where ω is the angular frequency of the given collective mode, the system of equations of motion reads
(26)-Eθϕsin⁡θδϕ-Eθθsin⁡θδθ-λ~δϕ=0,Eϕϕsin⁡θδϕ+Eϕθsin⁡θδθ-λ~δθ=0.
The linear and homogeneous system of equations expressed in (26) can be written as an eigenvalue problem
(27)C__v→=λ~v→,
with λ~=i(Ms/γ)ω the complex eigenvalues of the problem,
(28)C__=[-Eθϕsin⁡θ-Eθθsin⁡θEϕϕsin⁡θEϕθsin⁡θ],
being C__ a real, but not symmetric, matrix and v→=(δϕ,δθ)T.However, the system of motion equations (26) can be also recast as a complex generalized Hermitian eigenvalue problem
(29)A__v→=λB__v→,
where B__ is a Hessian matrix expressed by the second derivatives of the energy density at equilibrium. In particular
(30)B__=[EϕϕEϕθEθϕEθθ],A__=[0isin⁡θ-isin⁡θ0].
It should be noticed that the matrix A__ is Hermitian, whereas B__ is real and symmetric, so that all the corresponding eigenvalues λ=γ/Msω are real quantities.


### 2.2.2. HDMM for an Isolated System Composed byN Interacting Magnetic Momenta

Equation (24) can be generalized to the case of N interacting magnetic momenta (where each momentum is identified with a micromagnetic cell) taking the form
(31)E=E0+12∑n=1N∑l=1N[Eϕnϕlδϕnδϕl+2Eϕnθlδϕnδθl+Eθnθlδθnδθl].
In (31), the total energy density is given by E=Eext+Eexch+Edmg+Eani where the different contributions are expressed in (4), (6), (7), and (9), respectively. By substituting (31) into (23), we get
(32)δϕk.=-γMssin⁡θk∑l=1N[Eθkϕlδϕl+Eθkθlδθl],δθk.=γMssin⁡θk∑l=1N[Eϕkϕlδϕl+Eϕkθlδθl].
By introducing the time dependence in the form eiωt, the system of equations of motion is composed by the following 2Nlinear and homogeneous equations for k=1⋯N:
(33)∑l=1N(-Eθkϕlsin⁡θk)δϕl+∑l=1N(-Eθkθlsin⁡θk)δθl-λ~δϕk=0,∑l=1N(Eϕkϕlsin⁡θk)δϕl+∑l=1N(Eϕkθlsin⁡θk)δθl-λ~δθk=0.
The unknown factors δϕl,δθl represent the eigenvectors of the problem and are expressed by the small angular deviation from the equilibrium position of the azimuthal (ϕl) and polar (θl) angles in the lth micromagnetic cell. The system above has a solution only if the determinant is zero. By suitable exchanges of rows (columns), the linear and homogeneous system of equations expressed in (33) can be written as an eigenvalue problem in analogy with the case of a macrospin system
(34)C__v→=λ~v→.
In (34), v→ is the set of the unknown factors representing the eigenvectors of the problem that take the form
(35)v→=(δϕ1,δθ1,δϕ2,δθ2,…,δϕN,δθN)T.C__ is the matrix whose elements are expressed as
(36)C2k-1,2l-1=-Eθkϕlsin⁡θkC2k-1,2l=-Eθkθlsin⁡θkC2k,2l-1=Eϕkϕlsin⁡θkC2k,2l=Eϕkθlsin⁡θk}k=1⋯N,l=1⋯N.
Analogously to the case of a macrospin system, the matrix C__ is real but not symmetric, also in this case the eigenvalues (and not only the eigenvectors) are complex. This matrix can be indeed seen as composed by two submatrices 2×2 for each pair of values (k,l). In the diagonal submatrices (k=l) the following relation is verified:
(37)C2k-1,2k-1=-C2k,2k.
For the elements of two different submatrices (k,l) and (l,k), that are not diagonal (k≠l), the following symmetries hold:
(38)sin⁡θkC2k-1,2l-1=-sin⁡θlC2l,2k,sin⁡θkC2k-1,2l=sin⁡θlC2l-1,2k,sin⁡θkC2k,2l-1=sin⁡θlC2l,2k-1,sin⁡θkC2k,2l=-sin⁡θlC2l-1,2k-1.
As for the macrospin approximation, the equation of motion can be recast as a complex generalized Hermitian eigenvalue problem
(39)A__v→=λB__v→,
where B__ is a Hessian matrix expressed by the second derivatives of the energy density at equilibrium. B__ is given by
(40)B2k-1,2l-1=EϕkϕlB2k-1,2l=EϕkθlB2k,2l-1=EθkϕlB2k,2l=Eθkθl}k=1⋯N,l=1⋯N,
where the matrix B__ is, again, real and symmetric. Moreover, since the static magnetization corresponds to a minimum of the energy and the matrix B__ is its Hessian, the matrix B__ is also positive defined. Instead, the matrix A__ has the following form:
(41)A__=[0isin⁡θ100⋯-isin⁡θ1000⋯000isin⁡θ2⋯00-isin⁡θ20⋯⋯⋯].
The matrix A__ is Hermitian. This allows us to solve the system as a complex generalized Hermitian eigenvalue problem which admits only real eigenvalues. To further reduce the computational time, it is possible to evaluate only some eigenvalues and eigenvectors that are in a specific range.Once the eigenvectorsv→ are obtained, the dynamic magnetization δm→k in the kth micromagnetic cell expressed in Cartesian coordinates and in unit of MS is given by
(42)δm→k=(-sin⁡θksin⁡ϕkδϕk+cos⁡θkcos⁡ϕkδθk,sisin⁡θkcos⁡ϕkδϕk+cos⁡θksin⁡ϕkδθk,-sin⁡θkδθk).
For each solution of the eigenvalue problem, the collection of all δm→k defines the mode profile. It must be remarked that δm→k is a complex vector, because δθk,δϕk are, in general, complex.


## 2.2.1. HDMM for a Macrospin System

The equation of motion for the magnetic systems under investigation will be derived by following semiclassical approach. Also, the model will be first derived by considering the so-called macrospin approximation, where the material is thought as uniformly magnetized and represented by a single dipole momentum.As known from classical mechanics, in the presence of fixed constraints and conservative sources, the system Hamiltonianℋ coincides with the total mechanical energy E~, namely, ℋ≡E~=T-U, where T is the kinetic energy and U is the potential expressed as the opposite of the potential energy V. By defining the Lagrangian variables of the problem with qn, where n=1,2,… is the number of degrees of freedom corresponding to the dynamic variables, and the corresponding conjugate momenta with pn, the Hamilton equations in the 2n canonical variables (qn, pn) take the form [15]
(19)∂qn∂t=∂H∂pn,∂pn∂t=-∂H∂qn.
The direction of the magnetic dipole moment of thekth cell is given by (3). For the specific case, the dynamic variables referred to the k-cell are the small deviations from equilibrium of the azimuthal and polar angles given by
(20)q1=δϕ,q2=δθ,
where 1 (2) labels the first (second) variable.To determine the conjugate momenta, the expression of the angular momentum is needed. We recall the relation between the angular momentuml→Ω and the magnetic momentum μ→ by referring to a a rigid body with a fixed point Ω (see Figure 1), namely,
(21)l→Ω=1γμ→=vcMsγm→,
where γ is the gyromagnetic ratio and vc the volume of the magnetic moment to the case. Indeed, since q1 represents a rotation about thez-axis of the magnetic dipole, its conjugate momentum p1 corresponds to the z component of the variation of the angular momentum, namely,(22a)p1=δlz=vcMsγδmz=-vcMsγsin⁡θδθ,
where m→=(sinθcosϕ,sinθsinϕ,cosθ) is the unit magnetization vector, normalized to the saturation magnetization Ms. Following an analogous argument, the other momentum p2 can be determined. Indeed, q2 is a rotation of the dipole moment about an axis of unitary vector ϕ^′=-sinϕe^1+cosϕe^2 in the x-yplane at an angle ϕ^′=ϕ+π/2 from the x-axis. Therefore, p2 corresponds to the projection of the variation of the angular momentum along the ϕ^′ vector, namely,(22b)p2=δl→Ω⋅ϕ′^=vcMsγ(δm→⋅ϕ^′)=-vcMsγ(∂m→∂ϕδϕ+∂m→∂θδθ)⋅ϕ^′=vcMsγsin⁡θδϕ.By substituting (20), (22a), and (22b) into the first Hamilton equation (cf. (19)), the following system of equations is obtained:
(23)δϕ.=-γvcMssin⁡θHδθ,δθ.=γvcMssin⁡θHδϕ,
where the dot notation stands for the time derivative and, at the right-hand side, the first derivatives of the energy with respect to the mechanical variables appear. By introducing the energy density E=E~/V (keeping in mind that ℋ=E~) and expanding it in a power Taylor expansion around the equilibrium up to the second order, it yields
(24)E=E0+12[Eϕϕ(δϕ)2+2Eϕθδϕδθ+Eθθ(δθ)2],
where E0 is the constant zero-order term that is inessential, the first-order terms vanish at equilibrium, and Eαβ represent the second derivatives calculated at equilibrium (Eαβ=∂2E/∂δα∂δβ with α,β=ϕ,θ). By using (23) and (24), we obtain
(25)δϕ.=-γMssin⁡θ[Eθϕδϕ+Eθθδθ],δθ.=γMssin⁡θ[Eϕϕδϕ+Eϕθδθ].
By inserting the time dependence in the form eiωt, where ω is the angular frequency of the given collective mode, the system of equations of motion reads
(26)-Eθϕsin⁡θδϕ-Eθθsin⁡θδθ-λ~δϕ=0,Eϕϕsin⁡θδϕ+Eϕθsin⁡θδθ-λ~δθ=0.
The linear and homogeneous system of equations expressed in (26) can be written as an eigenvalue problem
(27)C__v→=λ~v→,
with λ~=i(Ms/γ)ω the complex eigenvalues of the problem,
(28)C__=[-Eθϕsin⁡θ-Eθθsin⁡θEϕϕsin⁡θEϕθsin⁡θ],
being C__ a real, but not symmetric, matrix and v→=(δϕ,δθ)T.However, the system of motion equations (26) can be also recast as a complex generalized Hermitian eigenvalue problem
(29)A__v→=λB__v→,
where B__ is a Hessian matrix expressed by the second derivatives of the energy density at equilibrium. In particular
(30)B__=[EϕϕEϕθEθϕEθθ],A__=[0isin⁡θ-isin⁡θ0].
It should be noticed that the matrix A__ is Hermitian, whereas B__ is real and symmetric, so that all the corresponding eigenvalues λ=γ/Msω are real quantities.


## 2.2.2. HDMM for an Isolated System Composed byN Interacting Magnetic Momenta

Equation (24) can be generalized to the case of N interacting magnetic momenta (where each momentum is identified with a micromagnetic cell) taking the form
(31)E=E0+12∑n=1N∑l=1N[Eϕnϕlδϕnδϕl+2Eϕnθlδϕnδθl+Eθnθlδθnδθl].
In (31), the total energy density is given by E=Eext+Eexch+Edmg+Eani where the different contributions are expressed in (4), (6), (7), and (9), respectively. By substituting (31) into (23), we get
(32)δϕk.=-γMssin⁡θk∑l=1N[Eθkϕlδϕl+Eθkθlδθl],δθk.=γMssin⁡θk∑l=1N[Eϕkϕlδϕl+Eϕkθlδθl].
By introducing the time dependence in the form eiωt, the system of equations of motion is composed by the following 2Nlinear and homogeneous equations for k=1⋯N:
(33)∑l=1N(-Eθkϕlsin⁡θk)δϕl+∑l=1N(-Eθkθlsin⁡θk)δθl-λ~δϕk=0,∑l=1N(Eϕkϕlsin⁡θk)δϕl+∑l=1N(Eϕkθlsin⁡θk)δθl-λ~δθk=0.
The unknown factors δϕl,δθl represent the eigenvectors of the problem and are expressed by the small angular deviation from the equilibrium position of the azimuthal (ϕl) and polar (θl) angles in the lth micromagnetic cell. The system above has a solution only if the determinant is zero. By suitable exchanges of rows (columns), the linear and homogeneous system of equations expressed in (33) can be written as an eigenvalue problem in analogy with the case of a macrospin system
(34)C__v→=λ~v→.
In (34), v→ is the set of the unknown factors representing the eigenvectors of the problem that take the form
(35)v→=(δϕ1,δθ1,δϕ2,δθ2,…,δϕN,δθN)T.C__ is the matrix whose elements are expressed as
(36)C2k-1,2l-1=-Eθkϕlsin⁡θkC2k-1,2l=-Eθkθlsin⁡θkC2k,2l-1=Eϕkϕlsin⁡θkC2k,2l=Eϕkθlsin⁡θk}k=1⋯N,l=1⋯N.
Analogously to the case of a macrospin system, the matrix C__ is real but not symmetric, also in this case the eigenvalues (and not only the eigenvectors) are complex. This matrix can be indeed seen as composed by two submatrices 2×2 for each pair of values (k,l). In the diagonal submatrices (k=l) the following relation is verified:
(37)C2k-1,2k-1=-C2k,2k.
For the elements of two different submatrices (k,l) and (l,k), that are not diagonal (k≠l), the following symmetries hold:
(38)sin⁡θkC2k-1,2l-1=-sin⁡θlC2l,2k,sin⁡θkC2k-1,2l=sin⁡θlC2l-1,2k,sin⁡θkC2k,2l-1=sin⁡θlC2l,2k-1,sin⁡θkC2k,2l=-sin⁡θlC2l-1,2k-1.
As for the macrospin approximation, the equation of motion can be recast as a complex generalized Hermitian eigenvalue problem
(39)A__v→=λB__v→,
where B__ is a Hessian matrix expressed by the second derivatives of the energy density at equilibrium. B__ is given by
(40)B2k-1,2l-1=EϕkϕlB2k-1,2l=EϕkθlB2k,2l-1=EθkϕlB2k,2l=Eθkθl}k=1⋯N,l=1⋯N,
where the matrix B__ is, again, real and symmetric. Moreover, since the static magnetization corresponds to a minimum of the energy and the matrix B__ is its Hessian, the matrix B__ is also positive defined. Instead, the matrix A__ has the following form:
(41)A__=[0isin⁡θ100⋯-isin⁡θ1000⋯000isin⁡θ2⋯00-isin⁡θ20⋯⋯⋯].
The matrix A__ is Hermitian. This allows us to solve the system as a complex generalized Hermitian eigenvalue problem which admits only real eigenvalues. To further reduce the computational time, it is possible to evaluate only some eigenvalues and eigenvectors that are in a specific range.Once the eigenvectorsv→ are obtained, the dynamic magnetization δm→k in the kth micromagnetic cell expressed in Cartesian coordinates and in unit of MS is given by
(42)δm→k=(-sin⁡θksin⁡ϕkδϕk+cos⁡θkcos⁡ϕkδθk,sisin⁡θkcos⁡ϕkδϕk+cos⁡θksin⁡ϕkδθk,-sin⁡θkδθk).
For each solution of the eigenvalue problem, the collection of all δm→k defines the mode profile. It must be remarked that δm→k is a complex vector, because δθk,δϕk are, in general, complex.


## 2.3. Equations of Motion for Interacting Magnetic Particles

Let us suppose to have a 2D periodic array of interacting nanodots characterized by the primitive vectorsa→1 and a→2; for example, for the specific case of a rectangular lattice their values are
(43)a→1=λxx^,a→2=λyy^,
where λx and λy represent the periodicity along x-axis and y-axis, respectively. The primitive vectors of the reciprocal lattice are b→1 and b→2(44)b→1=2π(a→2×a→1)×a→2(a→1×a→2)2,b→2=2π(a→1×a→2)×a→2(a→1×a→2)2.
For the special case of a rectangular lattice the vectors become
(45)b→1=2πλxx^,b→2=2πλyy^.
Due to the analogy with the Bloch wave (analogy, not equality, because the wave function has not a physical meaning, contrarily to the magnetization), it is possible to write the following periodicity rule valid for the dynamic magnetization:
(46)δm→(r→+R→)=eiK→⋅R→δm→(r→),
where R→ is a vector of the particle lattice given by
(47)R→=i1a→1+i2a→2i1,i2∈Z,i1,i2=-Ni2⋯Ni2-1,
and r→ can be confined into the first primitive cell centered in the first dot. The Bloch vector takes the following values:
(48)K→=n1N1b→1+n2N2b→2,ni∈Z,ni=-Ni2⋯Ni2-1.N1,N2∈ℕ indicate the number of primitive cells n in direction a→1 and a→2, respectively. In this scheme, both N1 and N2 are taken as even numbers. In order to confirm the hypothesis on the dynamic magnetization N1,N2 must be very large.If the magnetizations of different primitive cells and the different micromagnetic cells were independent, that is no periodicity rule were present, then one would have a dynamic system with variablesθkR→ and ϕkR→, where the k index changes inside the magnetic particle and R→ can assume the values indicated in (47). In this case, the linear and homogeneous system of 2N equations of motion is
(49)∑l,R→′(-EθkR→ϕlR→′sin⁡θkR→)δϕlR→′+∑l,R′→(-EθkR→θlR→′sin⁡θjkR→)δθlR→′-λδϕkR→=0,∑l,R→′(EϕkR→ϕlR→′sin⁡θkR→)δϕlR→′+∑l,R′→(EϕkR→θlR→′sin⁡θkR→)δθlR→′-λδθkR→=0,
where k=1⋯N and the sums over l and R→′ are on the same values. Instead, thanks to the Bloch condition expressed in (46), one can consider the equations only at R→=0; moreover, taking into account the same condition, the variables appearing for R→′≠0 can be replaced by using the same condition. Now, when rewriting the system, the index R→ is omitted when it has value equal to 0 or it is irrelevant. Owing to these considerations, the system given in (49) can be rewritten in the form [16]
(50)∑l=1N(-∑R→′EθkϕlR→′eiK→⋅R→′sin⁡θk)δϕl+∑l=1N(-∑R→′EθkθlR→′eiK→⋅R→′sin⁡θk)δθlδθl-λ~δϕk=0,∑l=1N(∑R→′EϕkϕlR→′eiK→⋅R→′sin⁡θk)δϕl+∑l=1N(∑R→′EϕkθlR→′eiK→⋅R→′sin⁡θk)δθlδθl-λ~δθk=0.
Equation (50) is similar to (33) by making the following replacement:
(51)Eαkβl⟶∑R→′eiK→⋅R→′Eαk0βlR→′,
and recalling that now the energy is referred to the whole system of particles. Like for the case of the isolated magnetic particle, also for the case of interacting magnetic particles the system of linear and homogeneous equations given in (50) can be written as an eigenvalue problem which in turn can be cast as a complex generalized Hermitian eigenvalue problem.The symmetry for the matrix elements that was valid for a single primitive cell now is not respected except forK→=0 or K→=G→/2 with G→ a translational reciprocal vector:
(52)∑R→′eiK→⋅R→′Eαk0→βlR→′=∑R→′eiK→⋅R→′EβlR→′αk0→=∑R→′eiK→⋅R→′Eβl0→αk-R→′=∑R→′e-iK→⋅R→″Eβl0→αkR→″≠∑R→′eiK→⋅R→′Eβl0→αkR→′,K→≠0→,K→≠G→2,
with R→″=-R→′. The primitive cell has at the centre a single dot that occupies only a part of it. The interdot exchange coupling is zero. Thanks to the last consideration and to the fact that derivatives of Zeeman, exchange, anisotropy, energy density are referred only to the cell of the first variable (αk0→) with αkR→=θkR→,ϕkR→ or at most to the nearest neighbour, all terms of the sum in (51) with R→′≠0 are zero. Hence, for these energy density terms, the equations are the same as those of the single particle case and the same occurs for their corresponding derivatives appearing in the equations of motion. The only energy density term that differs from the one obtained for the isolated element is the demagnetizing energy density. For a system of interacting nanoparticles, the demagnetizing energy density can be written as
(53)Edmg=12μ0∑R→,R′→,k,k′m→k(R→)N__(R→,R→′,k,k′)m→k′(R→′).
Due to its rather complex expression, it is useful to give the derivation also of the first derivative of Edmg like for the case of the isolated nanoelement.In order to calculate the first derivative, the properties of the demagnetizing tensor must be considered(54a)〈v|N|ω〉=〈ω|N+|v〉*=〈ω|N+|v〉=〈ω|Nt|v〉=〈ω|N|v〉,
since N is real and symmetric and(54b)N__(R→,R→′,k,k′)=N__(R→′-R→+r→k′-r→k)=N__(-R→′+R→-r→k′+r→k)=N__(R→′,R→,k′,k),
thanks to the inversion symmetry. Hence, the first derivative of the energy is(55)∂Edmg∂δαk(0→)μ0=μ0Ms22[∑R→′k′(R→′,k′)≠(0→,k)m→k′(R→′)⋅N__(R→′,0→,k′,k)∂m→k(0→)∂δαk(0→)μ0=μ0Ms22μ0+∑R→′k′(R→′,k′)≠(0→,k)∂m→k(0→)∂δαk(0→)⋅N__(0→,R→′,k′,k)m→k′(R→′)μ0=μ0Ms22μ0+∂∂δαk(0→)(m→k(0→)⋅N__(0→,0→,k,k)m→k(0→))∑R→′k′(R→′,k′)≠(0→,k)m→k′(R→′)⋅N__(R→′,0→,k′,k)∂m→k(0→)∂δαk(0→)]μ0=μ0Ms2[∑R→′k′(R→′,k′)≠(0→,k)m→k′(R→′)⋅N__(R→′,0→,k′,k)∂m→k(0→)∂δαk(0→)μ0=μ0Ms22μ0+m→k(0→)⋅N__(0→,0→,k,k)∂m→k(0→)∂δαk(0→)∑R→′k′(R→′,k′)≠(0→,k)m→k′(R→′)⋅N__(R→′,0→,k′,k)∂m→k(0→)∂δαk(0→)].
The second derivative takes the form
(56)∂2Edmg∂δαk(0→)∂δβk(R→)μ0={μ0MS2[∑R→′k′m→k′(R→′)⋅N__(R→′,0→,k′,k)μ0MS2μ0×∂2m→k(0→)∂δαk(0→)∂δβk(0→)μ0MS2μ0+∂m→k(0→)∂δβk(0→)μ0MS2μ0⋅N__(0→,0→,k,k)∂m→k(0→)∂δαk(0→)]R→=0→,l=kμ0Ms2∂m→l(0→)∂δβl(0→)⋅N__(R→,0→,l,k)∂m→k(0→)∂δαk(0→)(R→,l)≠(0→,k).
In the sum it is included the case in which (R→′,k′)=(0→,k) that is generated from the derivative of the second term in the expression of the first derivative given above. Now it is possible to calculate the terms that enter into the system of (50), starting with the one corresponding to l=k(57a)∑R→eiK→⋅R→Eαk0→βlR→=∑R→eiK→⋅R→Eαk0→βkR→∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2[∑R→′k′m→k′(0→)⋅N__(R→′,0→,k′,k)∑R→eiK→⋅R→Eαk0→βlR→×∂2m→k(0→)∂δαk(0→)∂δβk(0→)+∂m→k(0→)∂δβk(0→)∑R→eiK→⋅R→Eαk0→βlR→⋅N__(0→,0→,k,k)∂m→k(0→)∂δαk(0→)]∑R→eiK→⋅R→Eαk0→βlR→+μ0Ms2∑R→≠0→eiK→⋅R→∂m→k(0→)∂δβk(0→)∑R→eiK→⋅R→Eαk0→βlR→⋅N__(R→,0→,k,k)∂m→k(0→)∂δαk(0→)∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2∑R→[∑k′m→k′(0→)⋅N__(R→,0→,k′,k)∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2μ0Ms2×∂2m→k(0→)∂δαk(0→)∂δβk(0→)∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2μ0Ms2+eiK→⋅R→∂m→k(0→)∂δβk(0→)∑R→eiK→⋅R→Eαk0→βlR→=μ0Ms2μ0Ms2⋅N__(R→,0→,k,k)∂m→k(0→)∂δαk(0→)].
Equation (57a) is obtained by taking into account (46), by including the term corresponding to R→=0 into the sum performed over R→≠0 and by considering that both the static magnetization and the dynamic magnetization do not depend on R→. Indeed, second derivatives are calculated at equilibrium and the exponential eiK→·R→ appears on both the numerator and the denominator of the derivative.Whenl≠k, the term turns out to be
(57b)∑R→eiK→⋅R→Eαk0→βlR→μ0=μ0Ms2[∑R→eiK→⋅R→∂m→l(0→)∂δβl(0→)⋅N__(R→,0→,l,k)∂m→k(0→)∂δαk(0→)].Due to the properties of the demagnetizing tensor, the symmetry(58)∑R→′eiK→⋅R→′Eαk0→βlR→′=∑R→′eiK→⋅R→′Eαl0→βkR→′,
is fulfilled when l=k, but it is not fulfilled when l≠k.The formalism previously developed for interacting particles can be extended to a system of 2D antidots (ADs). In this case, it is necessary to add the exchange interaction between primitive cells [21]. In extended magnetic system like AD arrays, in addition to the usual nearest-neighbours exchange interaction between micromagnetic cells, the exchange contribution across the nearest-neighbours micromagnetic cells belonging to adjacent surface primitive cells must be taken into account. Hence, we recall the exchange energy density of (6)
(59)Eexch=A∑k∑n(1-m→k⋅m→n)akn2.
Here, the first sum runs over all the micromagnetic cells of the primitive cell and the second sum runs over the nearest neighbours of the kth micromagnetic cell. When the kth micromagnetic cell is on one of the edges (vertices) of the given primitive cell, the interaction with one (two) micromagnetic cell(s) belonging to the correct nearest primitive cell must be added.


## 3. LDMM

This section is devoted to the review of the LDMM approach through which we derive the generalized Lagrange equation in the presence of two dissipative effects arising from the “positive” intrinsic damping and the “negative” one induced by the current-driven spin-transfer torque [22]. As for the Hamiltonian approach, we limit our study to the dynamics taking place in the linear and autonomous regime.


### 3.1. Description of the Magnetic System and Equation of Motion

The magnetic systems in which such competing phenomena take place are generally referred to as spin-valve nanopillars. These are heterostructures composed by two ferromagnetic layers, having generally different thicknesses, separated by a nonmagnetic (metallic or insulating) spacer, which is used to decouple the exchange interactions between them. The thinner magnetic layer is generally referred to as “Free Layer” (FL), whereas the thicker one is called “Fixed Layer” or “Pinned Layer” (PL). By means of an external voltage source and metallic contacts applied at the top and bottom of the structure, a current flow traverses the structure along the normal-to-plane direction (see Figure2).Figure 2
A schematic of a spin-valve nanopillar device.In this kind of device, conservative effects arise from the previously mentioned classical micromagnetic contributions (exchange, demagnetizing, magnetocrystalline anisotropy, and Zeeman fields) together with the Ampere (or Oersted) field due to the current flow. This latter contribution, however, will be neglected for simplicity. In fact, the main goal of this section is to describe why and how a Lagrangian approach needs to be take into account when nonconservative forces act in the system. The inclusion of more sophisticated effects, such as the Ampere field, will become relevant when this approach will be further generalized to describe the dynamics occurring in the nonlinear regime of spin-wave generation.Let us now briefly recall the governing equation of motion. When no dissipative contributions are taken into account, a persistent precessional motion of the magnetization vectorM→ takes place. It is described by classical Landau-Lifshitz equation M→.=γ(H→eff×M→), where H→eff is the effective field which accounts for all the above-mentioned contributions.On the other hand, nonconservative contributions arise from the material intrinsic dissipation and the spin-transfer torque induced by the current flow [26]. The former accounts for the phenomenologically introduced intrinsic Gilbert dissipation [27], the phenomenon by which the precessional motion of M→, excited by a given stimulus, relaxes towards its equilibrium state. The relaxation rate is proportional to a scalar quantity, called Gilbert constant α. The torque exerted on the magnetization is generally represented by T→ID=(α/MS)(M→×M→.).Concerning the dissipative effects induced by the current flow, it has been extensively shown that this bias current can become spin-polarized in the direction of the magnetization vector of the thicker magnetic layer and can then transfer this induced spin angular momentum to the magnetization of the thinner magnetic layer. For a proper direction of the bias currentI, this spin-transfer mechanism creates a torque which opposes to that induced by the Gilbert damping, creating an effective negative damping. The corresponding spin-transfer torque, derived by Slonczewski [28], can be expressed as T→ST=(σI/MS)[M→×(M→×p→)], where the unit vector p→ defines the direction of the spin polarization (in turn defined by the magnetization vector of the PL) and the constant σ modulates the strength of the spin-torque effect. It is equal to σ=εg0μB/2eMSSd, where ε is the spin-torque efficiency (defined in [29]), g0 is the Landè factor, μB is the Bohr magneton, e is the absolute value of the electronic charge, d is the thickness of the magnetic layer, and S is the current-carrying area [30].It should be mentioned that, because of the larger value of both saturation magnetization and thickness, the PL is not substantially affected by any current-driven magnetization dynamics, so that it is generally treated as it were fixed (or pinned) along its equilibrium direction. On the contrary, the FL’s properties allow it to describe more easily several kinds of dynamics (e.g., switching [31], precession [32], domain-wall motion [33], gyrotropic motion of vortex state [34]), as it were “free” to move.Under these circumstances, the magnetization dynamics of the FL is governed by the Landau-Lifshitz-Gilbert-Slonczewski (LLGS) equation(60)M→.=γ(H→eff×M→)+αMS(M→×M→.)+σIMS[M→×(M→×p→)].
It has to be remarked that the equality (with opposite sign) of the two torques (T→ID=-T→ST), achievable by means of a proper intensity and sign of current (the positive one, I > 0, which corresponds to a current flow moving from the PL to the FL), which in turn implies the fully compensation of the two dissipation mechanisms, yields the system in an out-of-equilibrium zero-dissipation stationary state (a limit cycle, using the notation of dynamic systems) (see Figure 3). In such a regime, the excitation of microwave spin waves becomes physically conceivable.Figure 3
Schematic representation of undamped (in the absence of both Gilbert and current-induced damping) and out-of-equilibrium (in the presence of both Gilbert and current-induced damping) precessions. In this latter case, a limit cycle is described as well since the torques due to the intrinsic dissipation and the spin-transfer-induced one (for a proper intensity and direction of the current) balance each other.To derive the mathematical formulation of the LDMM, for the sake of simplicity, we write, first, the generalized Lagrange equation for the case of an isolated magnetic particle within the macrospin approximation. After that, we will generalize this approach for the case of an isolated particle composed byN interacting momenta.


#### 3.1.1. LDMM for a Macrospin System

We will preliminarily assume that the dynamics of the magnetization vector could be described through three degress of freedom(θ,ϕ,ψ), as shown in Figure 1. In such a framework the generalized Lagrange equations read
(61)ddt(∂L∂ϕ.)-∂L∂ϕ+∂IST∂ϕ.+∂IID∂ϕ.=0,ddt(∂L∂θ.)-∂L∂θ+∂IST∂θ.+∂IID∂θ.=0,ddt(∂L∂ψ.)-∂L∂ψ+∂IST∂ψ.+∂IID∂ψ.=0,
where L = T+ U represents the Lagrangian of the system given by the sum of the kinetic energy T and the potential U, whereas ℑST and ℑID are the dissipation functions related to the spin torque and the intrinsic damping, respectively.Taking advantage of the explicit formulations of the energy contributions given in Section2.1, we can rewrite the previous system by accounting for the relationship among the potential, the kinetic energy, and the conservative part E~=T-U of the total mechanical energy (which in the conservative limit coincides with the Hamiltonian of the system ℋ in the presence of fixed constraints). By substituting L with 2T-E~, we thus obtain
(62)2ddt(∂T∂ϕ.)-ddt(∂E~∂ϕ.)-2∂T∂ϕ+∂E~∂ϕ+∂IST∂ϕ.+∂IID∂ϕ.=0,2ddt(∂T∂θ.)-ddt(∂E~∂θ.)-2∂T∂θ+∂E~∂θ+∂IST∂θ.+∂IID∂θ.=0,2ddt(∂T∂ψ.)-ddt(∂E~∂ψ.)-2∂T∂ψ+∂E~∂ψ+∂IST∂ψ.+∂IID∂ψ.=0.
To solve the proposed problem, we need thus to explicitly find the expressions of the variables T, E~, ℑID and ℑST, shown in (62), in Lagrangian coordinates.The kinetic energyT associated to the precessional motion of the magnetization vector can be computed by referring to the case of a rigid body with a fixed point Ω (see Figure 1). In this case, the rotational kinetic energy density is expressed as
(63)T=12l→Ω⋅ω→,
where l→Ω is the angular momentum and ω→ is the angular velocity vector. Notice that in the following the kinetic energy T, the dissipation functions ℑID and ℑST, the work dw, and the energy losses dw/dt are referred to their corresponding quantities per unit volume. In addition, we recall the definition of the energy density E=E~/V given in Section 2.1. By virtue of the relationship between angular momentum l→Ω and magnetic dipole momentum μ→ (l→Ω=(1/γ)μ→), and considering that this latter is directed along e^3′ (see Figure 1), the previous scalar product involves, in turn, the only component of l→Ω along e^3′. It leads to
(64)T=12MSγ(ϕ.cos⁡θ+ψ.).
Equation (64) therefore represents the rotational kinetic energy T=T(θ,ϕ,ψ) expressed in Lagrangian coordinates with MS=μ/V. The conservative part E of the energy density accounts for all the standard micromagnetic contributions discussed previously. As shown in Section 2.1, all the contributions appearing in E only depend on the Lagrangian variables θ,ϕ, but not on their derivatives, namely, E=E(θ,ϕ).To derive the Lagrangian formulation of the dissipative contributions, let us start from the classical definition of the workdw carried out by a magnetic system subject to a nonconservative force. As it is known, such a force has to be derived from the gradient of a dissipation function ℑ and the rate of energy losses associated to a dissipative torque can be thus expressed as
(65)dwdt=±M→.⋅∂I∂M→.,
where the plus (minus) sign accounts for torques which act as “drain” (“source”) of energy and refers to the case ℑ = ℑID (ℑ = ℑST).So, by multiplying the LLGS equation (60) by H→eff and assuming that the energy losses rates are small compared to the conservative (precessional) part, namely, M→.≃-γμ0(M→×H→eff), we deduce the following expressions for the dissipative function densities
(66)IID=α2γM0M→.2,(67)IST=σJγM0p→⋅(M→×M→.).
It should be noticed that (66), which appears in the usual form of a Rayleigh-like dissipation function, is a positive-definite form, as expected for a power dissipated through a viscous friction mechanism, whatever the magnetization configuration (α, γ, and Ms are positive constants). On the other hand, the spin-torque dissipation function of (67) is a non-Rayleigh one and strongly depends, apart from the direction of the current flow J, on the relative magnetization configuration of the ferromagnetic layers (M→ and p→).Finally, taking into account the expression of the time-independent unit vectorp→=(sinΘPLcosΦPL,sinΘPLsinΦPL,cosΘPL), the explicit expressions of the dissipative functions densities are
(68)IID=αMs2γ(θ.2+ϕ.2sin2θ)=IID(θ,ϕ),IST=σJMsγ[f1(ϕ)θ.+f2(θ,ϕ)ϕ.]=IST(θ,ϕ),
where
(69)f1(ϕ)=sin⁡ΘPLsin⁡(ΦPL-ϕ),(70)f2(θ,ϕ)=sin2θcos⁡ΘPL-cos⁡θsin⁡θsin⁡ΘPLcos⁡(ΦPL-ϕ).
We proceed by evaluating first the third Lagrange equation
(71)2ddt(∂T∂ψ.)-ddt(∂E∂ψ.)-2∂T∂ψ+∂E∂ψ+∂IST∂ψ.+∂IID∂ψ.=0.
In (71), we notice that the kinetic energy is the only term dependent on the intrinsic rotation, and, in particular, on its velocity ψ·, so that the previous equation reduces to
(72)2ddt(∂T∂ψ.)=0,
which stands for a first integral of motion representing the conservation of the (only component of the) angular momentum
(73)2∂T∂ψ.=MSγ=constant≡l→Ω⋅e→3′.
Such a first integral also points out that our system can be described through only two degrees of freedom, as expected for the characterization of the dynamics of a vector having constant modulus which undergoes a precession (with a fixed point) onto a sphere. Taking into account this result, in the following we will use the parameters θ,ϕ.The first two Lagrange equations read(74)2ddt(∂T∂ϕ.)-ddt(∂E∂ϕ.)-2∂T∂ϕ+∂E∂ϕ+∂IST∂ϕ.+∂IID∂ϕ.=0,2ddt(∂T∂θ.)-ddt(∂E∂θ.)-2∂T∂θ+∂E∂θ+∂IST∂θ.+∂IID∂θ.=0,
which, substituting the corresponding expressions (the conservative part E of the total energy density will be discussed later on in the text), become
(75)-Msγsin⁡θθ.+∂E∂ϕ+σJMsγf2+αMsγϕ.sin2θ=0,Msγϕ.sin⁡θ+∂E∂θ+σJMsγf1+αMsγθ.=0.Since in the present LDMM approach we are interested in the characterization of the linear dynamics of these magnetic current-driven auto-oscillatory systems, we adopt the classical formalism of small oscillations. By following this procedure, we linearize the system of (75) by considering small perturbations (δθ, δϕ) around the equilibrium configuration of the FL (Θ,Φ), as done in (2), and expand the energy terms in Taylor series up to the second perturbative order.By means of such a formalism, the system given in (75) leads to the set of generalized Lagrange equations
(76)-Msγsin⁡Θδθ.+∂E∂δϕ+σJMsγMsγ×[χ(Θ,Φ)δθ+η(Θ,Φ)δϕ+ν(Θ,Φ)]Msγ-+αMsγδϕ.sin2Θ=0,Msγδϕ.sin⁡Θ+∂E∂δθ+σJMsγ[ζ(Θ,Φ)δϕ+β(Θ,Φ)]Msγ+αMsγδθ.=0,
where
(77)χ(Θ,Φ)=sin⁡2Θcos⁡ΘPL-cos⁡2Θsin⁡ΘPLcos⁡(ΦPL-Φ),η(Θ,Φ)=cos⁡Θsin⁡Θsin⁡ΘPLsin⁡(Φ-ΦPL),ν(Θ,Φ)=sin⁡Θ(sin⁡Θcos⁡ΘPL-cos⁡Θsin⁡ΘPLcos⁡(ΦPL-Φ)),ζ(Θ,Φ)=-sin⁡ΘPLcos⁡(ΦPL-Φ),β(Θ,Φ)=sin⁡ΘPLsin⁡(ΦPL-Φ).Furthermore, by developing also the energy densityE in Taylor series around the equilibrium state in analogy with what was done in (24) for HDMM, we notice that only the second derivatives appear in the equation. Indeed, the inessential constant term can be neglected and the first derivatives vanish at equilibrium. Moreover, considering that the new Lagrangian variables (δθ,δϕ) have to exhibit a time dependence proportional to exp(iωt) (being ω/2π the frequency of the spin-wave eigenmode), we end up with
(78)δϕ[Eϕϕ+σJMsγη]+δθ[Eϕθ+σJMsγχ]δϕ=δϕ[-iωαMsγsin2Θ]+δθ[iωMsγsin⁡Θ],δϕ[Eθϕ+σJMsγζ]+δθ[Eθθ]δϕ=δϕ[-iωMsγsin⁡Θ]+δθ[-iωαMsγ],
where the subscripts stand for partial derivative with respect to the indicated variables, whose explicit expressions can be found in the Section 2.1. The terms involving β(Θ,Φ) and ν(Θ,Φ) have been disregarded since they do not exhibit an explicit time dependence.By settingλ=γ/Msω and σ~=σMs/γ, (78) can be recast in the form of a complex generalized non-Hermitian eigenvalue problem
(79)A__v→=λB__v→,
with
(80)A__=[-iαsin2Θisin⁡Θ-isin⁡Θ-iα],B__=[Eϕϕ+σ~JηEϕθ+σ~JχEθϕ+σ~JζEθθ],v→=[δϕ,δθ]T.
As expected, if no damping and current are taken into account (α = 0, J = 0 A/m2), the system in (79) so obtained recovers exactly the fully conservative HDMM one, where the matrix A__ is Hermitian, whereas the matrix B__ is real and symmetric (see (40)), and all the corresponding eigenvalues are real.In the presence of dissipative effects due to damping and spin torque, the symmetry of the problem is strongly reduced and the corresponding eigenvalues will be, in general, complex quantities: the real part represents the mode frequency, whereas the imaginary part defines the mode decay rate. With this information in hands, it will be possible to establish, after a single iteration, the subset of the spin-wave normal modes which becomes unstable after the application of a spin-polarized current. In detail, for any applied current below the excitation threshold, the imaginary parts of all eigenvalues have to be positive, recalling the behavior of a damped oscillator. On the contrary, for current values above the excitation threshold, the imaginary parts of the activated normal modes switch to a negative value, giving rise to the instability mechanism which determines, in the time domain, the growth of the precession cone and, in turn, the change of the generated frequency.Finally, it should be mentioned that, unlike HDMM, LDMM cannot be recast in a nongeneralized eigenvalue problem form.


#### 3.1.2. LDMM for an Isolated System Composed byN Interacting Magnetic Momenta

Let us consider now the magnetic system of our micromagnetic calculations composed by a finite numberN of interacting magnetic momenta. Within our nanoscale numerical approach, each magnetic momentum is identified by means of a micromagnetic cell. The generalization of (78) to N cells (k=1,…,N) leads to a system of 2N scalar and linear equations as follows:
(81)λ{∑l=1N[Eϕkϕlδϕl+Eϕkθlδθl]+δθk[σ~Jχk]+δϕk[σ~Jηk]}λλ=δϕk[-iαsin2Θk]+δθk[isin⁡Θk],λ{∑l=1N[Eθkϕlδϕl+Eθkθlδθl]+δϕk[σ~Jζk]}λλ=δϕk[-isin⁡Θk]+δθk[-iα].
The system in (81) can be analogously formulated as a complex generalized non-Hermitian eigenvalue problem (in the same form as (79)), where
(82)A__=[-iαsin2Θ1isin⁡Θ100⋯-isin⁡Θ1-iα00⋯00-iαsin2Θ2isin⁡Θ1⋯00-isin⁡Θ2-iα⋯⋯⋯⋯⋯⋯],B__=[Eϕ1ϕ1+σ~Jη1Eϕ1θ1+σ~Jχ1Eϕ1ϕ2Eϕ1θ2⋯Eθ1ϕ1+σ~Jζ1Eθ1θ1Eθ1ϕ2Eθ1θ2⋯Eϕ2ϕ1Eϕ2θ1Eϕ2ϕ2+σ~Jη2Eϕ2θ2+σ~Jχ2⋯Eθ2ϕ1Eθ2θ1Eθ1ϕ1+σ~Jζ1Eθ2θ2⋯⋯⋯⋯⋯⋯],v→=[δϕ1,δθ1,δϕ2,δθ2,……δϕN,δθN]T.
It is also interesting to notice that both matrices A__ and B__ appearing in (80) and (82) admit a decomposition which allows to separate the conservative part from the nonconservative one
(83)A__=Ac__+Anc__B__=Bc__+Bnc__,
where
(84)Ac__=[0isin⁡Θ100⋯-isin⁡Θ1000⋯000isin⁡Θ1⋯00-isin⁡Θ20⋯⋯⋯⋯⋯⋯],Anc__=[-iαsin2Θ1000⋯0-iα00⋯00-iαsin2Θ20⋯000-iα⋯⋯⋯⋯⋯⋯],Bc__=[Eϕ1ϕ1Eϕ1θ1Eϕ1ϕ2Eϕ1θ2⋯Eθ1ϕ1Eθ1θ1Eθ1ϕ2Eθ1θ2⋯Eϕ2ϕ1Eϕ2θ1Eϕ2ϕ2Eϕ2θ2⋯Eθ2ϕ1Eθ2θ1Eθ1ϕ1Eθ2θ2⋯⋯⋯⋯⋯⋯],Bnc__=[σ~Jη1σ~Jχ100⋯σ~Jζ1000⋯00σ~Jη2σ~Jχ2⋯00σ~Jζ10⋯⋯⋯⋯⋯⋯].
The results of this approach have been successfully compared with those coming from another micromagnetic framework which integrates the LLGS equation in the time domain by using a finite-difference scheme [22, 23, 35]. In these works, we evaluated the accuracy of the LDMM approach in determining the excitation threshold and studied in detail the reorientational phase transition which takes place when the direction of the external magnetic field is varied.Finally, we would like to mention that the extension of the LDMM approach to model the more realistic (and attractive) nonlinear and nonautonomous dynamics [30, 36–41] is currently under study.


## 3.1. Description of the Magnetic System and Equation of Motion

The magnetic systems in which such competing phenomena take place are generally referred to as spin-valve nanopillars. These are heterostructures composed by two ferromagnetic layers, having generally different thicknesses, separated by a nonmagnetic (metallic or insulating) spacer, which is used to decouple the exchange interactions between them. The thinner magnetic layer is generally referred to as “Free Layer” (FL), whereas the thicker one is called “Fixed Layer” or “Pinned Layer” (PL). By means of an external voltage source and metallic contacts applied at the top and bottom of the structure, a current flow traverses the structure along the normal-to-plane direction (see Figure2).Figure 2
A schematic of a spin-valve nanopillar device.In this kind of device, conservative effects arise from the previously mentioned classical micromagnetic contributions (exchange, demagnetizing, magnetocrystalline anisotropy, and Zeeman fields) together with the Ampere (or Oersted) field due to the current flow. This latter contribution, however, will be neglected for simplicity. In fact, the main goal of this section is to describe why and how a Lagrangian approach needs to be take into account when nonconservative forces act in the system. The inclusion of more sophisticated effects, such as the Ampere field, will become relevant when this approach will be further generalized to describe the dynamics occurring in the nonlinear regime of spin-wave generation.Let us now briefly recall the governing equation of motion. When no dissipative contributions are taken into account, a persistent precessional motion of the magnetization vectorM→ takes place. It is described by classical Landau-Lifshitz equation M→.=γ(H→eff×M→), where H→eff is the effective field which accounts for all the above-mentioned contributions.On the other hand, nonconservative contributions arise from the material intrinsic dissipation and the spin-transfer torque induced by the current flow [26]. The former accounts for the phenomenologically introduced intrinsic Gilbert dissipation [27], the phenomenon by which the precessional motion of M→, excited by a given stimulus, relaxes towards its equilibrium state. The relaxation rate is proportional to a scalar quantity, called Gilbert constant α. The torque exerted on the magnetization is generally represented by T→ID=(α/MS)(M→×M→.).Concerning the dissipative effects induced by the current flow, it has been extensively shown that this bias current can become spin-polarized in the direction of the magnetization vector of the thicker magnetic layer and can then transfer this induced spin angular momentum to the magnetization of the thinner magnetic layer. For a proper direction of the bias currentI, this spin-transfer mechanism creates a torque which opposes to that induced by the Gilbert damping, creating an effective negative damping. The corresponding spin-transfer torque, derived by Slonczewski [28], can be expressed as T→ST=(σI/MS)[M→×(M→×p→)], where the unit vector p→ defines the direction of the spin polarization (in turn defined by the magnetization vector of the PL) and the constant σ modulates the strength of the spin-torque effect. It is equal to σ=εg0μB/2eMSSd, where ε is the spin-torque efficiency (defined in [29]), g0 is the Landè factor, μB is the Bohr magneton, e is the absolute value of the electronic charge, d is the thickness of the magnetic layer, and S is the current-carrying area [30].It should be mentioned that, because of the larger value of both saturation magnetization and thickness, the PL is not substantially affected by any current-driven magnetization dynamics, so that it is generally treated as it were fixed (or pinned) along its equilibrium direction. On the contrary, the FL’s properties allow it to describe more easily several kinds of dynamics (e.g., switching [31], precession [32], domain-wall motion [33], gyrotropic motion of vortex state [34]), as it were “free” to move.Under these circumstances, the magnetization dynamics of the FL is governed by the Landau-Lifshitz-Gilbert-Slonczewski (LLGS) equation(60)M→.=γ(H→eff×M→)+αMS(M→×M→.)+σIMS[M→×(M→×p→)].
It has to be remarked that the equality (with opposite sign) of the two torques (T→ID=-T→ST), achievable by means of a proper intensity and sign of current (the positive one, I > 0, which corresponds to a current flow moving from the PL to the FL), which in turn implies the fully compensation of the two dissipation mechanisms, yields the system in an out-of-equilibrium zero-dissipation stationary state (a limit cycle, using the notation of dynamic systems) (see Figure 3). In such a regime, the excitation of microwave spin waves becomes physically conceivable.Figure 3
Schematic representation of undamped (in the absence of both Gilbert and current-induced damping) and out-of-equilibrium (in the presence of both Gilbert and current-induced damping) precessions. In this latter case, a limit cycle is described as well since the torques due to the intrinsic dissipation and the spin-transfer-induced one (for a proper intensity and direction of the current) balance each other.To derive the mathematical formulation of the LDMM, for the sake of simplicity, we write, first, the generalized Lagrange equation for the case of an isolated magnetic particle within the macrospin approximation. After that, we will generalize this approach for the case of an isolated particle composed byN interacting momenta.


### 3.1.1. LDMM for a Macrospin System

We will preliminarily assume that the dynamics of the magnetization vector could be described through three degress of freedom(θ,ϕ,ψ), as shown in Figure 1. In such a framework the generalized Lagrange equations read
(61)ddt(∂L∂ϕ.)-∂L∂ϕ+∂IST∂ϕ.+∂IID∂ϕ.=0,ddt(∂L∂θ.)-∂L∂θ+∂IST∂θ.+∂IID∂θ.=0,ddt(∂L∂ψ.)-∂L∂ψ+∂IST∂ψ.+∂IID∂ψ.=0,
where L = T+ U represents the Lagrangian of the system given by the sum of the kinetic energy T and the potential U, whereas ℑST and ℑID are the dissipation functions related to the spin torque and the intrinsic damping, respectively.Taking advantage of the explicit formulations of the energy contributions given in Section2.1, we can rewrite the previous system by accounting for the relationship among the potential, the kinetic energy, and the conservative part E~=T-U of the total mechanical energy (which in the conservative limit coincides with the Hamiltonian of the system ℋ in the presence of fixed constraints). By substituting L with 2T-E~, we thus obtain
(62)2ddt(∂T∂ϕ.)-ddt(∂E~∂ϕ.)-2∂T∂ϕ+∂E~∂ϕ+∂IST∂ϕ.+∂IID∂ϕ.=0,2ddt(∂T∂θ.)-ddt(∂E~∂θ.)-2∂T∂θ+∂E~∂θ+∂IST∂θ.+∂IID∂θ.=0,2ddt(∂T∂ψ.)-ddt(∂E~∂ψ.)-2∂T∂ψ+∂E~∂ψ+∂IST∂ψ.+∂IID∂ψ.=0.
To solve the proposed problem, we need thus to explicitly find the expressions of the variables T, E~, ℑID and ℑST, shown in (62), in Lagrangian coordinates.The kinetic energyT associated to the precessional motion of the magnetization vector can be computed by referring to the case of a rigid body with a fixed point Ω (see Figure 1). In this case, the rotational kinetic energy density is expressed as
(63)T=12l→Ω⋅ω→,
where l→Ω is the angular momentum and ω→ is the angular velocity vector. Notice that in the following the kinetic energy T, the dissipation functions ℑID and ℑST, the work dw, and the energy losses dw/dt are referred to their corresponding quantities per unit volume. In addition, we recall the definition of the energy density E=E~/V given in Section 2.1. By virtue of the relationship between angular momentum l→Ω and magnetic dipole momentum μ→ (l→Ω=(1/γ)μ→), and considering that this latter is directed along e^3′ (see Figure 1), the previous scalar product involves, in turn, the only component of l→Ω along e^3′. It leads to
(64)T=12MSγ(ϕ.cos⁡θ+ψ.).
Equation (64) therefore represents the rotational kinetic energy T=T(θ,ϕ,ψ) expressed in Lagrangian coordinates with MS=μ/V. The conservative part E of the energy density accounts for all the standard micromagnetic contributions discussed previously. As shown in Section 2.1, all the contributions appearing in E only depend on the Lagrangian variables θ,ϕ, but not on their derivatives, namely, E=E(θ,ϕ).To derive the Lagrangian formulation of the dissipative contributions, let us start from the classical definition of the workdw carried out by a magnetic system subject to a nonconservative force. As it is known, such a force has to be derived from the gradient of a dissipation function ℑ and the rate of energy losses associated to a dissipative torque can be thus expressed as
(65)dwdt=±M→.⋅∂I∂M→.,
where the plus (minus) sign accounts for torques which act as “drain” (“source”) of energy and refers to the case ℑ = ℑID (ℑ = ℑST).So, by multiplying the LLGS equation (60) by H→eff and assuming that the energy losses rates are small compared to the conservative (precessional) part, namely, M→.≃-γμ0(M→×H→eff), we deduce the following expressions for the dissipative function densities
(66)IID=α2γM0M→.2,(67)IST=σJγM0p→⋅(M→×M→.).
It should be noticed that (66), which appears in the usual form of a Rayleigh-like dissipation function, is a positive-definite form, as expected for a power dissipated through a viscous friction mechanism, whatever the magnetization configuration (α, γ, and Ms are positive constants). On the other hand, the spin-torque dissipation function of (67) is a non-Rayleigh one and strongly depends, apart from the direction of the current flow J, on the relative magnetization configuration of the ferromagnetic layers (M→ and p→).Finally, taking into account the expression of the time-independent unit vectorp→=(sinΘPLcosΦPL,sinΘPLsinΦPL,cosΘPL), the explicit expressions of the dissipative functions densities are
(68)IID=αMs2γ(θ.2+ϕ.2sin2θ)=IID(θ,ϕ),IST=σJMsγ[f1(ϕ)θ.+f2(θ,ϕ)ϕ.]=IST(θ,ϕ),
where
(69)f1(ϕ)=sin⁡ΘPLsin⁡(ΦPL-ϕ),(70)f2(θ,ϕ)=sin2θcos⁡ΘPL-cos⁡θsin⁡θsin⁡ΘPLcos⁡(ΦPL-ϕ).
We proceed by evaluating first the third Lagrange equation
(71)2ddt(∂T∂ψ.)-ddt(∂E∂ψ.)-2∂T∂ψ+∂E∂ψ+∂IST∂ψ.+∂IID∂ψ.=0.
In (71), we notice that the kinetic energy is the only term dependent on the intrinsic rotation, and, in particular, on its velocity ψ·, so that the previous equation reduces to
(72)2ddt(∂T∂ψ.)=0,
which stands for a first integral of motion representing the conservation of the (only component of the) angular momentum
(73)2∂T∂ψ.=MSγ=constant≡l→Ω⋅e→3′.
Such a first integral also points out that our system can be described through only two degrees of freedom, as expected for the characterization of the dynamics of a vector having constant modulus which undergoes a precession (with a fixed point) onto a sphere. Taking into account this result, in the following we will use the parameters θ,ϕ.The first two Lagrange equations read(74)2ddt(∂T∂ϕ.)-ddt(∂E∂ϕ.)-2∂T∂ϕ+∂E∂ϕ+∂IST∂ϕ.+∂IID∂ϕ.=0,2ddt(∂T∂θ.)-ddt(∂E∂θ.)-2∂T∂θ+∂E∂θ+∂IST∂θ.+∂IID∂θ.=0,
which, substituting the corresponding expressions (the conservative part E of the total energy density will be discussed later on in the text), become
(75)-Msγsin⁡θθ.+∂E∂ϕ+σJMsγf2+αMsγϕ.sin2θ=0,Msγϕ.sin⁡θ+∂E∂θ+σJMsγf1+αMsγθ.=0.Since in the present LDMM approach we are interested in the characterization of the linear dynamics of these magnetic current-driven auto-oscillatory systems, we adopt the classical formalism of small oscillations. By following this procedure, we linearize the system of (75) by considering small perturbations (δθ, δϕ) around the equilibrium configuration of the FL (Θ,Φ), as done in (2), and expand the energy terms in Taylor series up to the second perturbative order.By means of such a formalism, the system given in (75) leads to the set of generalized Lagrange equations
(76)-Msγsin⁡Θδθ.+∂E∂δϕ+σJMsγMsγ×[χ(Θ,Φ)δθ+η(Θ,Φ)δϕ+ν(Θ,Φ)]Msγ-+αMsγδϕ.sin2Θ=0,Msγδϕ.sin⁡Θ+∂E∂δθ+σJMsγ[ζ(Θ,Φ)δϕ+β(Θ,Φ)]Msγ+αMsγδθ.=0,
where
(77)χ(Θ,Φ)=sin⁡2Θcos⁡ΘPL-cos⁡2Θsin⁡ΘPLcos⁡(ΦPL-Φ),η(Θ,Φ)=cos⁡Θsin⁡Θsin⁡ΘPLsin⁡(Φ-ΦPL),ν(Θ,Φ)=sin⁡Θ(sin⁡Θcos⁡ΘPL-cos⁡Θsin⁡ΘPLcos⁡(ΦPL-Φ)),ζ(Θ,Φ)=-sin⁡ΘPLcos⁡(ΦPL-Φ),β(Θ,Φ)=sin⁡ΘPLsin⁡(ΦPL-Φ).Furthermore, by developing also the energy densityE in Taylor series around the equilibrium state in analogy with what was done in (24) for HDMM, we notice that only the second derivatives appear in the equation. Indeed, the inessential constant term can be neglected and the first derivatives vanish at equilibrium. Moreover, considering that the new Lagrangian variables (δθ,δϕ) have to exhibit a time dependence proportional to exp(iωt) (being ω/2π the frequency of the spin-wave eigenmode), we end up with
(78)δϕ[Eϕϕ+σJMsγη]+δθ[Eϕθ+σJMsγχ]δϕ=δϕ[-iωαMsγsin2Θ]+δθ[iωMsγsin⁡Θ],δϕ[Eθϕ+σJMsγζ]+δθ[Eθθ]δϕ=δϕ[-iωMsγsin⁡Θ]+δθ[-iωαMsγ],
where the subscripts stand for partial derivative with respect to the indicated variables, whose explicit expressions can be found in the Section 2.1. The terms involving β(Θ,Φ) and ν(Θ,Φ) have been disregarded since they do not exhibit an explicit time dependence.By settingλ=γ/Msω and σ~=σMs/γ, (78) can be recast in the form of a complex generalized non-Hermitian eigenvalue problem
(79)A__v→=λB__v→,
with
(80)A__=[-iαsin2Θisin⁡Θ-isin⁡Θ-iα],B__=[Eϕϕ+σ~JηEϕθ+σ~JχEθϕ+σ~JζEθθ],v→=[δϕ,δθ]T.
As expected, if no damping and current are taken into account (α = 0, J = 0 A/m2), the system in (79) so obtained recovers exactly the fully conservative HDMM one, where the matrix A__ is Hermitian, whereas the matrix B__ is real and symmetric (see (40)), and all the corresponding eigenvalues are real.In the presence of dissipative effects due to damping and spin torque, the symmetry of the problem is strongly reduced and the corresponding eigenvalues will be, in general, complex quantities: the real part represents the mode frequency, whereas the imaginary part defines the mode decay rate. With this information in hands, it will be possible to establish, after a single iteration, the subset of the spin-wave normal modes which becomes unstable after the application of a spin-polarized current. In detail, for any applied current below the excitation threshold, the imaginary parts of all eigenvalues have to be positive, recalling the behavior of a damped oscillator. On the contrary, for current values above the excitation threshold, the imaginary parts of the activated normal modes switch to a negative value, giving rise to the instability mechanism which determines, in the time domain, the growth of the precession cone and, in turn, the change of the generated frequency.Finally, it should be mentioned that, unlike HDMM, LDMM cannot be recast in a nongeneralized eigenvalue problem form.


### 3.1.2. LDMM for an Isolated System Composed byN Interacting Magnetic Momenta

Let us consider now the magnetic system of our micromagnetic calculations composed by a finite numberN of interacting magnetic momenta. Within our nanoscale numerical approach, each magnetic momentum is identified by means of a micromagnetic cell. The generalization of (78) to N cells (k=1,…,N) leads to a system of 2N scalar and linear equations as follows:
(81)λ{∑l=1N[Eϕkϕlδϕl+Eϕkθlδθl]+δθk[σ~Jχk]+δϕk[σ~Jηk]}λλ=δϕk[-iαsin2Θk]+δθk[isin⁡Θk],λ{∑l=1N[Eθkϕlδϕl+Eθkθlδθl]+δϕk[σ~Jζk]}λλ=δϕk[-isin⁡Θk]+δθk[-iα].
The system in (81) can be analogously formulated as a complex generalized non-Hermitian eigenvalue problem (in the same form as (79)), where
(82)A__=[-iαsin2Θ1isin⁡Θ100⋯-isin⁡Θ1-iα00⋯00-iαsin2Θ2isin⁡Θ1⋯00-isin⁡Θ2-iα⋯⋯⋯⋯⋯⋯],B__=[Eϕ1ϕ1+σ~Jη1Eϕ1θ1+σ~Jχ1Eϕ1ϕ2Eϕ1θ2⋯Eθ1ϕ1+σ~Jζ1Eθ1θ1Eθ1ϕ2Eθ1θ2⋯Eϕ2ϕ1Eϕ2θ1Eϕ2ϕ2+σ~Jη2Eϕ2θ2+σ~Jχ2⋯Eθ2ϕ1Eθ2θ1Eθ1ϕ1+σ~Jζ1Eθ2θ2⋯⋯⋯⋯⋯⋯],v→=[δϕ1,δθ1,δϕ2,δθ2,……δϕN,δθN]T.
It is also interesting to notice that both matrices A__ and B__ appearing in (80) and (82) admit a decomposition which allows to separate the conservative part from the nonconservative one
(83)A__=Ac__+Anc__B__=Bc__+Bnc__,
where
(84)Ac__=[0isin⁡Θ100⋯-isin⁡Θ1000⋯000isin⁡Θ1⋯00-isin⁡Θ20⋯⋯⋯⋯⋯⋯],Anc__=[-iαsin2Θ1000⋯0-iα00⋯00-iαsin2Θ20⋯000-iα⋯⋯⋯⋯⋯⋯],Bc__=[Eϕ1ϕ1Eϕ1θ1Eϕ1ϕ2Eϕ1θ2⋯Eθ1ϕ1Eθ1θ1Eθ1ϕ2Eθ1θ2⋯Eϕ2ϕ1Eϕ2θ1Eϕ2ϕ2Eϕ2θ2⋯Eθ2ϕ1Eθ2θ1Eθ1ϕ1Eθ2θ2⋯⋯⋯⋯⋯⋯],Bnc__=[σ~Jη1σ~Jχ100⋯σ~Jζ1000⋯00σ~Jη2σ~Jχ2⋯00σ~Jζ10⋯⋯⋯⋯⋯⋯].
The results of this approach have been successfully compared with those coming from another micromagnetic framework which integrates the LLGS equation in the time domain by using a finite-difference scheme [22, 23, 35]. In these works, we evaluated the accuracy of the LDMM approach in determining the excitation threshold and studied in detail the reorientational phase transition which takes place when the direction of the external magnetic field is varied.Finally, we would like to mention that the extension of the LDMM approach to model the more realistic (and attractive) nonlinear and nonautonomous dynamics [30, 36–41] is currently under study.


## 3.1.1. LDMM for a Macrospin System

We will preliminarily assume that the dynamics of the magnetization vector could be described through three degress of freedom(θ,ϕ,ψ), as shown in Figure 1. In such a framework the generalized Lagrange equations read
(61)ddt(∂L∂ϕ.)-∂L∂ϕ+∂IST∂ϕ.+∂IID∂ϕ.=0,ddt(∂L∂θ.)-∂L∂θ+∂IST∂θ.+∂IID∂θ.=0,ddt(∂L∂ψ.)-∂L∂ψ+∂IST∂ψ.+∂IID∂ψ.=0,
where L = T+ U represents the Lagrangian of the system given by the sum of the kinetic energy T and the potential U, whereas ℑST and ℑID are the dissipation functions related to the spin torque and the intrinsic damping, respectively.Taking advantage of the explicit formulations of the energy contributions given in Section2.1, we can rewrite the previous system by accounting for the relationship among the potential, the kinetic energy, and the conservative part E~=T-U of the total mechanical energy (which in the conservative limit coincides with the Hamiltonian of the system ℋ in the presence of fixed constraints). By substituting L with 2T-E~, we thus obtain
(62)2ddt(∂T∂ϕ.)-ddt(∂E~∂ϕ.)-2∂T∂ϕ+∂E~∂ϕ+∂IST∂ϕ.+∂IID∂ϕ.=0,2ddt(∂T∂θ.)-ddt(∂E~∂θ.)-2∂T∂θ+∂E~∂θ+∂IST∂θ.+∂IID∂θ.=0,2ddt(∂T∂ψ.)-ddt(∂E~∂ψ.)-2∂T∂ψ+∂E~∂ψ+∂IST∂ψ.+∂IID∂ψ.=0.
To solve the proposed problem, we need thus to explicitly find the expressions of the variables T, E~, ℑID and ℑST, shown in (62), in Lagrangian coordinates.The kinetic energyT associated to the precessional motion of the magnetization vector can be computed by referring to the case of a rigid body with a fixed point Ω (see Figure 1). In this case, the rotational kinetic energy density is expressed as
(63)T=12l→Ω⋅ω→,
where l→Ω is the angular momentum and ω→ is the angular velocity vector. Notice that in the following the kinetic energy T, the dissipation functions ℑID and ℑST, the work dw, and the energy losses dw/dt are referred to their corresponding quantities per unit volume. In addition, we recall the definition of the energy density E=E~/V given in Section 2.1. By virtue of the relationship between angular momentum l→Ω and magnetic dipole momentum μ→ (l→Ω=(1/γ)μ→), and considering that this latter is directed along e^3′ (see Figure 1), the previous scalar product involves, in turn, the only component of l→Ω along e^3′. It leads to
(64)T=12MSγ(ϕ.cos⁡θ+ψ.).
Equation (64) therefore represents the rotational kinetic energy T=T(θ,ϕ,ψ) expressed in Lagrangian coordinates with MS=μ/V. The conservative part E of the energy density accounts for all the standard micromagnetic contributions discussed previously. As shown in Section 2.1, all the contributions appearing in E only depend on the Lagrangian variables θ,ϕ, but not on their derivatives, namely, E=E(θ,ϕ).To derive the Lagrangian formulation of the dissipative contributions, let us start from the classical definition of the workdw carried out by a magnetic system subject to a nonconservative force. As it is known, such a force has to be derived from the gradient of a dissipation function ℑ and the rate of energy losses associated to a dissipative torque can be thus expressed as
(65)dwdt=±M→.⋅∂I∂M→.,
where the plus (minus) sign accounts for torques which act as “drain” (“source”) of energy and refers to the case ℑ = ℑID (ℑ = ℑST).So, by multiplying the LLGS equation (60) by H→eff and assuming that the energy losses rates are small compared to the conservative (precessional) part, namely, M→.≃-γμ0(M→×H→eff), we deduce the following expressions for the dissipative function densities
(66)IID=α2γM0M→.2,(67)IST=σJγM0p→⋅(M→×M→.).
It should be noticed that (66), which appears in the usual form of a Rayleigh-like dissipation function, is a positive-definite form, as expected for a power dissipated through a viscous friction mechanism, whatever the magnetization configuration (α, γ, and Ms are positive constants). On the other hand, the spin-torque dissipation function of (67) is a non-Rayleigh one and strongly depends, apart from the direction of the current flow J, on the relative magnetization configuration of the ferromagnetic layers (M→ and p→).Finally, taking into account the expression of the time-independent unit vectorp→=(sinΘPLcosΦPL,sinΘPLsinΦPL,cosΘPL), the explicit expressions of the dissipative functions densities are
(68)IID=αMs2γ(θ.2+ϕ.2sin2θ)=IID(θ,ϕ),IST=σJMsγ[f1(ϕ)θ.+f2(θ,ϕ)ϕ.]=IST(θ,ϕ),
where
(69)f1(ϕ)=sin⁡ΘPLsin⁡(ΦPL-ϕ),(70)f2(θ,ϕ)=sin2θcos⁡ΘPL-cos⁡θsin⁡θsin⁡ΘPLcos⁡(ΦPL-ϕ).
We proceed by evaluating first the third Lagrange equation
(71)2ddt(∂T∂ψ.)-ddt(∂E∂ψ.)-2∂T∂ψ+∂E∂ψ+∂IST∂ψ.+∂IID∂ψ.=0.
In (71), we notice that the kinetic energy is the only term dependent on the intrinsic rotation, and, in particular, on its velocity ψ·, so that the previous equation reduces to
(72)2ddt(∂T∂ψ.)=0,
which stands for a first integral of motion representing the conservation of the (only component of the) angular momentum
(73)2∂T∂ψ.=MSγ=constant≡l→Ω⋅e→3′.
Such a first integral also points out that our system can be described through only two degrees of freedom, as expected for the characterization of the dynamics of a vector having constant modulus which undergoes a precession (with a fixed point) onto a sphere. Taking into account this result, in the following we will use the parameters θ,ϕ.The first two Lagrange equations read(74)2ddt(∂T∂ϕ.)-ddt(∂E∂ϕ.)-2∂T∂ϕ+∂E∂ϕ+∂IST∂ϕ.+∂IID∂ϕ.=0,2ddt(∂T∂θ.)-ddt(∂E∂θ.)-2∂T∂θ+∂E∂θ+∂IST∂θ.+∂IID∂θ.=0,
which, substituting the corresponding expressions (the conservative part E of the total energy density will be discussed later on in the text), become
(75)-Msγsin⁡θθ.+∂E∂ϕ+σJMsγf2+αMsγϕ.sin2θ=0,Msγϕ.sin⁡θ+∂E∂θ+σJMsγf1+αMsγθ.=0.Since in the present LDMM approach we are interested in the characterization of the linear dynamics of these magnetic current-driven auto-oscillatory systems, we adopt the classical formalism of small oscillations. By following this procedure, we linearize the system of (75) by considering small perturbations (δθ, δϕ) around the equilibrium configuration of the FL (Θ,Φ), as done in (2), and expand the energy terms in Taylor series up to the second perturbative order.By means of such a formalism, the system given in (75) leads to the set of generalized Lagrange equations
(76)-Msγsin⁡Θδθ.+∂E∂δϕ+σJMsγMsγ×[χ(Θ,Φ)δθ+η(Θ,Φ)δϕ+ν(Θ,Φ)]Msγ-+αMsγδϕ.sin2Θ=0,Msγδϕ.sin⁡Θ+∂E∂δθ+σJMsγ[ζ(Θ,Φ)δϕ+β(Θ,Φ)]Msγ+αMsγδθ.=0,
where
(77)χ(Θ,Φ)=sin⁡2Θcos⁡ΘPL-cos⁡2Θsin⁡ΘPLcos⁡(ΦPL-Φ),η(Θ,Φ)=cos⁡Θsin⁡Θsin⁡ΘPLsin⁡(Φ-ΦPL),ν(Θ,Φ)=sin⁡Θ(sin⁡Θcos⁡ΘPL-cos⁡Θsin⁡ΘPLcos⁡(ΦPL-Φ)),ζ(Θ,Φ)=-sin⁡ΘPLcos⁡(ΦPL-Φ),β(Θ,Φ)=sin⁡ΘPLsin⁡(ΦPL-Φ).Furthermore, by developing also the energy densityE in Taylor series around the equilibrium state in analogy with what was done in (24) for HDMM, we notice that only the second derivatives appear in the equation. Indeed, the inessential constant term can be neglected and the first derivatives vanish at equilibrium. Moreover, considering that the new Lagrangian variables (δθ,δϕ) have to exhibit a time dependence proportional to exp(iωt) (being ω/2π the frequency of the spin-wave eigenmode), we end up with
(78)δϕ[Eϕϕ+σJMsγη]+δθ[Eϕθ+σJMsγχ]δϕ=δϕ[-iωαMsγsin2Θ]+δθ[iωMsγsin⁡Θ],δϕ[Eθϕ+σJMsγζ]+δθ[Eθθ]δϕ=δϕ[-iωMsγsin⁡Θ]+δθ[-iωαMsγ],
where the subscripts stand for partial derivative with respect to the indicated variables, whose explicit expressions can be found in the Section 2.1. The terms involving β(Θ,Φ) and ν(Θ,Φ) have been disregarded since they do not exhibit an explicit time dependence.By settingλ=γ/Msω and σ~=σMs/γ, (78) can be recast in the form of a complex generalized non-Hermitian eigenvalue problem
(79)A__v→=λB__v→,
with
(80)A__=[-iαsin2Θisin⁡Θ-isin⁡Θ-iα],B__=[Eϕϕ+σ~JηEϕθ+σ~JχEθϕ+σ~JζEθθ],v→=[δϕ,δθ]T.
As expected, if no damping and current are taken into account (α = 0, J = 0 A/m2), the system in (79) so obtained recovers exactly the fully conservative HDMM one, where the matrix A__ is Hermitian, whereas the matrix B__ is real and symmetric (see (40)), and all the corresponding eigenvalues are real.In the presence of dissipative effects due to damping and spin torque, the symmetry of the problem is strongly reduced and the corresponding eigenvalues will be, in general, complex quantities: the real part represents the mode frequency, whereas the imaginary part defines the mode decay rate. With this information in hands, it will be possible to establish, after a single iteration, the subset of the spin-wave normal modes which becomes unstable after the application of a spin-polarized current. In detail, for any applied current below the excitation threshold, the imaginary parts of all eigenvalues have to be positive, recalling the behavior of a damped oscillator. On the contrary, for current values above the excitation threshold, the imaginary parts of the activated normal modes switch to a negative value, giving rise to the instability mechanism which determines, in the time domain, the growth of the precession cone and, in turn, the change of the generated frequency.Finally, it should be mentioned that, unlike HDMM, LDMM cannot be recast in a nongeneralized eigenvalue problem form.


## 3.1.2. LDMM for an Isolated System Composed byN Interacting Magnetic Momenta

Let us consider now the magnetic system of our micromagnetic calculations composed by a finite numberN of interacting magnetic momenta. Within our nanoscale numerical approach, each magnetic momentum is identified by means of a micromagnetic cell. The generalization of (78) to N cells (k=1,…,N) leads to a system of 2N scalar and linear equations as follows:
(81)λ{∑l=1N[Eϕkϕlδϕl+Eϕkθlδθl]+δθk[σ~Jχk]+δϕk[σ~Jηk]}λλ=δϕk[-iαsin2Θk]+δθk[isin⁡Θk],λ{∑l=1N[Eθkϕlδϕl+Eθkθlδθl]+δϕk[σ~Jζk]}λλ=δϕk[-isin⁡Θk]+δθk[-iα].
The system in (81) can be analogously formulated as a complex generalized non-Hermitian eigenvalue problem (in the same form as (79)), where
(82)A__=[-iαsin2Θ1isin⁡Θ100⋯-isin⁡Θ1-iα00⋯00-iαsin2Θ2isin⁡Θ1⋯00-isin⁡Θ2-iα⋯⋯⋯⋯⋯⋯],B__=[Eϕ1ϕ1+σ~Jη1Eϕ1θ1+σ~Jχ1Eϕ1ϕ2Eϕ1θ2⋯Eθ1ϕ1+σ~Jζ1Eθ1θ1Eθ1ϕ2Eθ1θ2⋯Eϕ2ϕ1Eϕ2θ1Eϕ2ϕ2+σ~Jη2Eϕ2θ2+σ~Jχ2⋯Eθ2ϕ1Eθ2θ1Eθ1ϕ1+σ~Jζ1Eθ2θ2⋯⋯⋯⋯⋯⋯],v→=[δϕ1,δθ1,δϕ2,δθ2,……δϕN,δθN]T.
It is also interesting to notice that both matrices A__ and B__ appearing in (80) and (82) admit a decomposition which allows to separate the conservative part from the nonconservative one
(83)A__=Ac__+Anc__B__=Bc__+Bnc__,
where
(84)Ac__=[0isin⁡Θ100⋯-isin⁡Θ1000⋯000isin⁡Θ1⋯00-isin⁡Θ20⋯⋯⋯⋯⋯⋯],Anc__=[-iαsin2Θ1000⋯0-iα00⋯00-iαsin2Θ20⋯000-iα⋯⋯⋯⋯⋯⋯],Bc__=[Eϕ1ϕ1Eϕ1θ1Eϕ1ϕ2Eϕ1θ2⋯Eθ1ϕ1Eθ1θ1Eθ1ϕ2Eθ1θ2⋯Eϕ2ϕ1Eϕ2θ1Eϕ2ϕ2Eϕ2θ2⋯Eθ2ϕ1Eθ2θ1Eθ1ϕ1Eθ2θ2⋯⋯⋯⋯⋯⋯],Bnc__=[σ~Jη1σ~Jχ100⋯σ~Jζ1000⋯00σ~Jη2σ~Jχ2⋯00σ~Jζ10⋯⋯⋯⋯⋯⋯].
The results of this approach have been successfully compared with those coming from another micromagnetic framework which integrates the LLGS equation in the time domain by using a finite-difference scheme [22, 23, 35]. In these works, we evaluated the accuracy of the LDMM approach in determining the excitation threshold and studied in detail the reorientational phase transition which takes place when the direction of the external magnetic field is varied.Finally, we would like to mention that the extension of the LDMM approach to model the more realistic (and attractive) nonlinear and nonautonomous dynamics [30, 36–41] is currently under study.



---

*Source: 765709-2012-08-13.xml*