# On the Formulation of the Exchange Field in the Landau-Lifshitz Equation for Spin-Wave Calculation in Magnonic Crystals

Authors: M. Krawczyk; M. L. Sokolovskyy; J. W. Klos; S. Mamica
Journal: Advances in Condensed Matter Physics
 (2012)

---

## Abstract

The calculation of the magnonic spectra using the plane-wave method has limitations, the origin of which lies in the formulation of the effective magnetic field term in the equation of motion (the Landau-Lifshitz equation) for composite media. According to ideas of the plane-wave method the system dynamics is described in terms of plane waves (a superposition of a number of plane waves), which are continuous functions and propagate throughout the medium. Since in magnonic crystals the sought-for superposition of plane waves represents the dynamic magnetization, the magnetic boundary conditions on the interfaces between constituent materials should be inherent in the Landau-Lifshitz equations. In this paper we present the derivation of the two expressions for the exchange field known from the literature. We start from the Heisenberg model and use a linear approximation and take into account the spacial dependence of saturation magnetization and exchange constant present in magnetic composites. We discuss the magnetic boundary conditions included in the presented formulations of the exchange field and elucidate their effect on spin-wave modes and their spectra in one- and two-dimensional planar magnonic crystals from plane-wave calculations.


---

## Body

## 1. Introduction

For the first time the exchange effects were discovered independently by W. Heisenberg and P.A.M. Dirac in 1926. They proposed the energy operator (Hamiltonian) for the exchange interaction between two particles with spinsS1 and S2 in the following form:
(1)Hex=-2J12S1⋅S2,
where J12 is the exchange integral. Because the exchange interaction is the fundamental one for magnetic materials than it is crucial for the calculations of the spin dynamics. An equation commonly used to describe the magnetization dynamics is the following Landau-Lifshitz (LL) equation:
(2)∂M(r,t)∂t=-γμ0[M(r,t)×(Hex(r,t)+⋯)]+⋯,
where Hex is the exchange field acting on the magnetization vector M. This equation is a macroscopic where all terms are in a form of continuous functions of a position vector r.The derivation of an exchange field in a uniform ferromagnetic material from the microscopic Heisenberg Hamiltonian (1) can be found in many textbooks, for example, [1–5]. We will follow ideas presented in these books but for composite materials, that is, when the structure consists of two or more constituent ferromagnetics being in direct contact.On the interface between two ferromagnetic materials, the boundary conditions (BCs) on dynamical component of the magnetization vector should be imposed. Such boundary conditions were proposed by Hoffman, then developed, and investigated by other authors [6–14]. From the LL equation together with the set of BCs, the spin-wave (SW) dispersion and profiles can be calculated. In many papers the calculation of the SW spectra in composite magnetic materials is based on the solution of the LL equation defined for a uniform material and then matched at the boundaries [7, 12, 13, 15–19]. In this paper we are interested in other method used to calculate the dispersion relation of the spin waves in magnetic composites with periodic distribution of constituent materials, and it is the plane-wave method (PWM). This method is widely used in calculations of the frequency spectra of an electromagnetic, elastic or electron waves propagating in a photonic crystal, phononic crystal or semiconductor periodic heterostructures, respectively. In this method, described in details in Section 3, boundary conditions at interfaces between constituent materials should be inherently included into the equation of motion, that is, by properly defined exchange field.Composites with a periodic arrangement of two (or more) different materials are extensively studied from many years. In the past, structures with periodicity in one dimension were investigated, these are multilayered structures which found many applications, for example, as a Bragg mirror or in GMR devices [20, 21]. Starting from 1987 and the discovery of photonic band gaps in photonic crystals [22, 23], the research was rapidly extended to other composites with periodicity in two and three dimensions. Among them there are phononic crystals, plasmonic crystals and also magnonic crystals (MCs) [24–30]. MCs can be regarded as magnetic analog of the photonic crystal, which uses the spin waves, instead of electromagnetic waves, to carry the information. MCs constitute one of the main building blocks of magnonic—promising direction of research focused on practical applications of spin waves [29–36]. For the development of magnonics, the computational methods have to be developed as well. The calculation of SW’s dispersion in MCs can be performed with different methods, for example, with micromagnetic simulations or dynamical matrix method but because of their complexity the computations are very time consuming [37–42]. It is important to develop other analytical and semianalytical methods, like a PWM which even though approximate will allow for efficient calculation of the dispersion of SWs in MCs with big insight into physical processes.Semiconductor periodic heterostructures (SHs) allow to tailor the electron and heat transport in nanoscale [43–45]. SHs are often described in effective mass approximation with the use of the envelope function instead of single electron functions. Effective mass equations derived from the Schrödinger equations are not unique, and many possible definitions of the kinetic energy operator were proposed. An extensive discussion about proper definition of the kinetic energy operator for the electron envelop function in SHs with position-dependent effective mass can be found in literature [46–58]. The calculation of the SW spectra in continuous model, that is, from the LL equation (2), undergoes similar difficulties, as we will show in this paper. However, this topic related to the calculations of the spin-wave dynamics in MCs is weakly presented in literature. In this paper we would like to fill this gap with a detailed consideration of different forms of the exchange field and then look at their consequences in the SW spectra of a MC calculated using the PWM.In this study we will show in details the derivation from the microscopic model different forms of the exchange field used for SW calculations in MCs. Then we will analyze differences in SW spectra in one- (1D) and two-dimensional (2D) thin films of MCs calculated with PWM for three different expressions of the exchange field. We will discuss the boundary condition implemented in each formulation. The paper consists of five sections. In Section2 we show the derivation of two forms of the exchange field from the Heisenberg Hamiltonian in linear approximation for magnetic composites with pointing at surface terms neglected. Then in Section 3 we introduce the PWM method and derive a final algebraic eigenvalue equations for different definitions of the exchange field. In Section 4 we present the results of the PWM calculations of SW spectra for these different forms of the exchange field for MCs. We will consider 1D and 2D MC. The paper finished with Section 5 where conclusions of our investigation are drawn.


## 2. Expression of the Exchange Field inInhomogeneous Media

We split the derivation of the expression of the exchange field in inhomogeneous materials into two steps. First, in Section2.1 we obtain the formula for the exchange energy density from the microscopic Heisenberg Hamiltonian. Here the crucial step is a transformation from the discrete model to the continuous one. In the second step, Section 2.2 the formula for the exchange field will be derived from the exchange energy density. In this step a linear approximation will be introduced, and the space dependence of magnetic material parameters will be considered.


### 2.1. Exchange Energy Functional

We start our calculations from the Heisenberg Hamiltonianℋl which defines exchange energy of the spin Sl on the lattice point l as follows:
(3)Hl=-2∑m∈(n.n.)‍JlmSl⋅Sm,
where Sm is the total spin vector on lattice point m, and the summation is performed over all nearest neighbors (n.n.) of a lth spin. Jlm is an exchange integral between the spins located at l and m. When we introduce normalized unit vector αl for the spin vector Sl(4)αl=Sl|Sl|,
then (3) will read
(5)Hl=-2|Sl|∑m∈(n.n.)‍Jlm|Sm|αl⋅αm.
According to the definition in (4), αl·αm=cosφ, where φ is an angle between spin vectors on lattice points l and m (see Figure 1).Figure 1
The discrete lattice of spins. The angle between neighboring spins:Sl and Sm is φ. We assume that φ and the spin deviation from the z axis are small.Let us assume that the angleφ between the nearest spin vectors is small and moreover that the spin vectors are continuous and smooth functions of a position vector r, that is, α=α(r) [2]. Formally we can do it through averaging S over the unit cell; that is, we introduce magnetization vector M(r) as follows:
(6)M(rl)=N¯μBgSl,whereN¯=NV
defines a number of spins (N) in the unit cell volume (V). μB is Bohr magneton, and g is a g factor (for free electrons g≈2). According to these definitions
(7)αl=Sl|Sl|≡M(rl)M(rl)
and we can write α(r) as a continuous function of the position vector as follows:
(8)α(r)=M(r)M(r).Having the continuous function of a position vector in hand, we can expand a unit vectorα in a lattice point m (α(rm)≡αm≡α′) in the Taylor series:
(9)α′=αl+∑i‍(∂xiαl)dxi+12∑i,j‍(∂xi∂xjαl)dxidxj+⋯,
where xi=x,y or z, and dxi is the distance between nearest spins along xi axis. ∂xi is an abbreviation of the partial derivative with respect to the Cartesian xi component. After limiting expansion up to the quadratic terms (it means that we assumed a small variation of α(r) in space) we can substitute this into (5) as follows:
(10)Hl=-2|Sl|∑m∈(n.n.)‍Jlm|Sm|-2|Sl|∑m∈(n.n.)‍Jlm|Sm|∑i‍αl⋅(∂xiαl)dxi-|Sl|∑m∈(n.n.)‍Jlm|Sm|∑i,j‍αl⋅(∂xi∂xjαl)dxidxj.
The second term on the right side is equal to zero because α is a unit vector, and the only possibility to change it is a rotation: unit vectors αl fulfill obvious relation: αl·αl=1. Differentiate this equation with respect to xi results in the following:
(11)(∂xiαl)⋅αl=0.
It means that ∂xiα is zero or is orthogonal to the vector α (i.e., αl·(∂xiαl)≡0, see (11)). The Hamiltonian can be now rewritten as follows:
(12)Hl=-2|Sl|∑m∈(n.n.)‍Jlm|Sm|-|Sl|∑m∈(n.n.)‍Jlm|Sm|∑i,j‍αl⋅(∂xi∂xjαl)dxidxj.For a homogeneous material the length of the spins is preserved,|Sl|=|Sm| in each lattice point, and the exchange integral is constant, that is, Jlm=J for each n.n. l and m (l≠m), (this assumption is valid for inhomogeneous material when two atomic planes (for n.n. exchange interactions) at the interfaces are removed from consideration.) With this homogenization, we obtain from (12) the following expression:
(13)Hl=-2Z|S|2J-|S|2J∑m∈(n.n.)‍∑i,j‍αl⋅(∂xi∂xjαl)dxidxj,
where Z is a number of the nearest neighbors. For crystals with cubic crystallographic structures (i.e., for simple cubic (sc), body-centered cubic (bcc), and face-centered cubic (fcc) lattice types) the distance between n.n. is equal along all directions, |dxi|=a. Using this property, the summation over nearest neighbors ∑m∈(n.n.)‍(⋯)dxidxj for i≠j is equal to 0 for each lattice type. So, we can obtain the following Hamiltonian:
(14)Hl=-2Z|S|2J-2|S|2Ja2∑i‍αl⋅(∂xi∂xiαl).Equation (11) can be again differentiate with respect to xj, and the result is
(15)αl⋅(∂xi∂xjαl)=-(∂xiαl)⋅(∂xjαl).
With this equality we can rewrite (14) for the exchange energy into the following form:
(16)Hl=-2Z|S|2J+2|S|2Ja2∑i(∂xiαl)2.To define energy density,ℰex as a continuous function of the position vector, we have to sum over all spins in the unit cell and divide it by the volume of this unit cell. In that way we obtain density of the following exchange energy:
(17)Eex=λM2+A∑i‍(∂xiα(r))2,
where
(18)A=2nJS2a
and n=1,2, or 4 for sc, bcc, or fcc lattice, respectively; [2], α(r) is defined in (8) and
(19)λ=-2ZJN¯μB2g2.To calculate the total exchange energy,Eex, stored in a magnetic material, we have to integrate density of the energy (17) over the volume of the material [3] as follows:
(20)Eex=∫V‍Eexd3r=∫V‍λM2d3r+∫V‍A∑i‍(∂xiα(r))2d3r=∫V‍λM2d3r+∫V‍A[(∂xMM)2+(∂yMM)2+(∂zMM)2]d3r.The SW can be regarded as coherent precession of the magnetization vector around its equilibrium direction. Based on this observation most SW calculation are performed in linear approximation. This approximation was already used once in our paper, it is in (10). Now we will use it again to simplify the expression (20) for the total exchange energy.


#### 2.1.1.  Exchange Energy in Form I

In linear approximation the magnetization vector(21)M(r)=Mx(r)x^+My(r)y^+Mz(r)z^
can be separated into two parts: a static and dynamic components. We assume that the magnetization component along the direction of the applied magnetic field, in our case it is the z-axis, is constant in time (but can be still position dependent), and its value is close to the length of the total magnetization vector as follows:
(22)M0=M0(r)≡Mz≈|M|≈const(t).
The time-depending components of the magnetization vector; Mx and My will be denoted by mx and my, respectively. We will define the dynamic magnetization vector as m=(mx,my) being a two-dimensional vector in the plane perpendicular to the direction of the saturation magnetization.The exchange energy, (20), with the help of approximation (22) can be rewritten in the following form:
(23)Eex=∫V‍λM02d3r+∫V‍A(∇mM0)2d3r.
This consists of the formula for the exchange energy, which we will call as Form I.


#### 2.1.2. Exchange Energy in Form II

In the following we will make further assumptions to obtain another expression for the exchange field. We can write that(24)(∂xmM0)2=((∂xm)M0-(∂xM0)mM02)2=(∂xm)2M02+(∂xM0)2m2-M0(∂xM0)(∂xm2)M04
were in the last component of the nominator we have used 2m·∂xm=∂x(m2). The same calculations can be applied to other components of the ∇ operator in the exchange energy functional (23). The following expression for the exchange energy can be obtained:
(25)Eex=∫V‍λM02d3r+∫V‍AM02(∇m)2d3r+∫V‍AM04(∇M0)2m2d3r-∫V‍AM03(∇M0)⋅(∇m2)d3r.
In MCs the saturation magnetization is a function of the position vector with a step increase at interfaces. For bicomponent MCs (i.e., consisting of two ferromagnetic materials: A and B), M0(r) can be defined with the help of the characteristic function S(r):
(26)M0(r)=(M0,A-M0,B)S(r)+M0,B,
were
(27)S(r)={1forrinmaterialA,0forrinmaterialB,
and M0,A andM0,B are saturation magnetizations in materials A and B, respectively. We can see that in last two terms in (25) there are derivatives of M0 with respect to the position which according to (26) are derivatives from the step function, that is, ∇M0=(M0,A-M0,B)δ(r-rinterface), where δ is the Dirac delta function, and rinterface is a position vector which define the interface. It means that these two terms are connected with the exchange energy contributed only at interfaces and which are related to the jump of the saturation magnetization value (in PWM calculations we will assume parallel magnetizations in both materials). It can be shown that these two terms result in internal magnetic field components localized on interfaces and that these components introduce singularities in the equation of motion. To avoid these singularities we neglect these two terms. (It can be shown by direct calculation of functional derivatives according to (29), that these terms introduce non-Hermitian (or non-anti-Hermitian) elements into equation of motion, that is, (m/M0)∇((2A/μ0M02)∇M0), and again its physical interpretation is questionable.)To summarize, we have derived two different formulas for the exchange energy in linear approximation, which are equivalent in the case of homogeneous material. These are(28)FormI:Eex=∫V‍λM02d3r+∫V‍A(∇mM0)2d3r,FormII:Eex=∫V‍λM02d3r+∫V‍AM02(∇m)2d3r.
In this derivation we have neglected the interface anisotropy terms [9, 11]. These effects, which can be present in real materials, have a microscopic origin and are limited to the very thin area around interfaces (one or two atomic planes). In continuous effective models such effects can be included by proper effective boundary conditions imposed on dynamical component of the magnetization vector.


### 2.2. Exchange Field

Exchange field can be derived from the exchange energy functionals (28) as a first variational derivative with respect to the magnetization vector [4, 59] as follows:
(29)Hex(r)=-1μ0δEexδM=-1μ0[δEexδmx,δEexδmy,δEexδM0].
This equation is written in SI units. Those variational derivatives can be calculated from Euler formula [4] as follows:
(30)δEexδmx=∂η∂mx-∇(∂η∂(∇mx))=∂η∂mx-∑i‍∂∂xi(∂η∂(∂ximx)),
where Eex=∫V‍ηd3r and η=λM02(r)+A(r)(∇m(r)/M0(r))2 or η=λM02(r)+(A(r)M02(r))(∇m(r))2 for Form I and Form II of the exchange energy, respectively (see (28)). We will perform those calculations independently for each form of the exchange energy defined in (28). During calculations we will take into account the inhomogeneity in the material, that is, spacial dependence of material parameters: A(r) and M0(r).


#### 2.2.1. Form I of the Exchange Field

For Form I of the exchange energy functionalη=λM02(r)+A(r)(∇(m(r)/M0(r)))2 and we will calculate functional derivative directly from (30). First we calculate the derivative with respect of m as follows:
(31)∂η∂mx=(∂A(r)∇(m(r)/M0(r)))2∂mx=A(r)M0(r)2∇(m(r)M0(r))⋅∇(M0(r)).
This term includes the derivative of saturation magnetization, which is a step function on the interface between two magnetic materials—(26). This part of the magnetic field is localized purely at interfaces similarly as was found already in (25). This term will introduce singularity into equation of motions, and it will be neglected.After this assumption the exchange magnetic field in the Form I will be obtained solely from the second term in (30), that is,
(32)μ0Hex=∇∂η∂∇m(r)=∇∂A(r)(∇(m(r)/M0(r)))2∂∇m(r)=∇[A(r)2(∇m(r)M0(r))⋅(∂∇(m(r)/M0(r))∂∇m(r))]=∇[2A(r)M0(r)]∇m(r)M0(r).


#### 2.2.2. Form II of the Exchange Field

For exchange energy written in the Form II as defined in (28), η=λM02(r)+(A(r)/M02(r))(∇m(r))2. We can calculate functional derivatives according to (30) and write Hex in a compact form without any further approximations. In this case ∂η/∂mx=0 because the first term in η is independent on ∂ximx. So the functional derivative of Eex with respect to mx take the following form:
(33)δEexδmx=-∑i‍∂∂xi(∂(A/M02)(∇m)2∂(∂ximx))=-∑i‍∂∂xi[AM02∂A∂(∂ximx)]=-∂∂x[2AM02(∂xmx)]-∂∂y[2AM02(∂ymx)]-∂∂z[2AM02(∂zmx)],
where 𝔄 denotes [(∂xmx)2+(∂xmy)2+(∂ymx)2+(∂ymy)2+(∂zmx)2+(∂zmy)2]. The same procedure can be repeated for the y component of the magnetization. Finally the exchange field can be written in the following form, that is, Form II:
(34)μ0Hex(r)=∇[2A(r)M02(r)]∇m(r).


#### 2.2.3. Summary of the Exchange Field Forms

We have shown the derivation of two different expressions for the exchange field in nonuniform ferromagnetic materials in linear approximation. These are(35)FormI:Hex(r)=∇lex,I(r)∇m~(r),wherem~(r)≡m(r)M0(r),lex,I=2A(r)μ0M0(r);FormII:Hex(r)=∇lex,II(r)∇m(r),wherelex,II=2A(r)μ0M02(r).We can also add to this list the exchange field in Form III, which is derived directly from the exchange energy functional (28) (independent of the form which will be used) under assumption of the homogeneous material, that is, when the space dependence of A and M0 is not taken into account during calculations of a functional derivative (30)
(36)FormIII:Hex(r)=lex,I(r)M0(r)∇2m(r).From the parameters introduced just above:lex,I and lex,II, only the second one (i.e., from Form II) has an additional physical meaning; that is, its square root defines the exchange length [3]. (In Form III the coefficient in the exchange field is the same as in Form II, but M0 was excluded from new parameter due to simplification in the latter calculus, see (42).)It is worth to note at this moment that differential operators in the definition of the exchange field, (35)-(36), work on dynamical components of the magnetization vector, m(r) in Forms II and III, while on its normalized function in Form I, m~(r). This will make different equations and should be kept in mind during the interpretation of the eigenvectors found in PWM and boundary conditions implemented in the exchange field definitions.Those three different formulas for the exchange field will be investigated for the calculation of the magnonic band structure in thin plates of 1D and 2D Mcs with PWM.


## 2.1. Exchange Energy Functional

We start our calculations from the Heisenberg Hamiltonianℋl which defines exchange energy of the spin Sl on the lattice point l as follows:
(3)Hl=-2∑m∈(n.n.)‍JlmSl⋅Sm,
where Sm is the total spin vector on lattice point m, and the summation is performed over all nearest neighbors (n.n.) of a lth spin. Jlm is an exchange integral between the spins located at l and m. When we introduce normalized unit vector αl for the spin vector Sl(4)αl=Sl|Sl|,
then (3) will read
(5)Hl=-2|Sl|∑m∈(n.n.)‍Jlm|Sm|αl⋅αm.
According to the definition in (4), αl·αm=cosφ, where φ is an angle between spin vectors on lattice points l and m (see Figure 1).Figure 1
The discrete lattice of spins. The angle between neighboring spins:Sl and Sm is φ. We assume that φ and the spin deviation from the z axis are small.Let us assume that the angleφ between the nearest spin vectors is small and moreover that the spin vectors are continuous and smooth functions of a position vector r, that is, α=α(r) [2]. Formally we can do it through averaging S over the unit cell; that is, we introduce magnetization vector M(r) as follows:
(6)M(rl)=N¯μBgSl,whereN¯=NV
defines a number of spins (N) in the unit cell volume (V). μB is Bohr magneton, and g is a g factor (for free electrons g≈2). According to these definitions
(7)αl=Sl|Sl|≡M(rl)M(rl)
and we can write α(r) as a continuous function of the position vector as follows:
(8)α(r)=M(r)M(r).Having the continuous function of a position vector in hand, we can expand a unit vectorα in a lattice point m (α(rm)≡αm≡α′) in the Taylor series:
(9)α′=αl+∑i‍(∂xiαl)dxi+12∑i,j‍(∂xi∂xjαl)dxidxj+⋯,
where xi=x,y or z, and dxi is the distance between nearest spins along xi axis. ∂xi is an abbreviation of the partial derivative with respect to the Cartesian xi component. After limiting expansion up to the quadratic terms (it means that we assumed a small variation of α(r) in space) we can substitute this into (5) as follows:
(10)Hl=-2|Sl|∑m∈(n.n.)‍Jlm|Sm|-2|Sl|∑m∈(n.n.)‍Jlm|Sm|∑i‍αl⋅(∂xiαl)dxi-|Sl|∑m∈(n.n.)‍Jlm|Sm|∑i,j‍αl⋅(∂xi∂xjαl)dxidxj.
The second term on the right side is equal to zero because α is a unit vector, and the only possibility to change it is a rotation: unit vectors αl fulfill obvious relation: αl·αl=1. Differentiate this equation with respect to xi results in the following:
(11)(∂xiαl)⋅αl=0.
It means that ∂xiα is zero or is orthogonal to the vector α (i.e., αl·(∂xiαl)≡0, see (11)). The Hamiltonian can be now rewritten as follows:
(12)Hl=-2|Sl|∑m∈(n.n.)‍Jlm|Sm|-|Sl|∑m∈(n.n.)‍Jlm|Sm|∑i,j‍αl⋅(∂xi∂xjαl)dxidxj.For a homogeneous material the length of the spins is preserved,|Sl|=|Sm| in each lattice point, and the exchange integral is constant, that is, Jlm=J for each n.n. l and m (l≠m), (this assumption is valid for inhomogeneous material when two atomic planes (for n.n. exchange interactions) at the interfaces are removed from consideration.) With this homogenization, we obtain from (12) the following expression:
(13)Hl=-2Z|S|2J-|S|2J∑m∈(n.n.)‍∑i,j‍αl⋅(∂xi∂xjαl)dxidxj,
where Z is a number of the nearest neighbors. For crystals with cubic crystallographic structures (i.e., for simple cubic (sc), body-centered cubic (bcc), and face-centered cubic (fcc) lattice types) the distance between n.n. is equal along all directions, |dxi|=a. Using this property, the summation over nearest neighbors ∑m∈(n.n.)‍(⋯)dxidxj for i≠j is equal to 0 for each lattice type. So, we can obtain the following Hamiltonian:
(14)Hl=-2Z|S|2J-2|S|2Ja2∑i‍αl⋅(∂xi∂xiαl).Equation (11) can be again differentiate with respect to xj, and the result is
(15)αl⋅(∂xi∂xjαl)=-(∂xiαl)⋅(∂xjαl).
With this equality we can rewrite (14) for the exchange energy into the following form:
(16)Hl=-2Z|S|2J+2|S|2Ja2∑i(∂xiαl)2.To define energy density,ℰex as a continuous function of the position vector, we have to sum over all spins in the unit cell and divide it by the volume of this unit cell. In that way we obtain density of the following exchange energy:
(17)Eex=λM2+A∑i‍(∂xiα(r))2,
where
(18)A=2nJS2a
and n=1,2, or 4 for sc, bcc, or fcc lattice, respectively; [2], α(r) is defined in (8) and
(19)λ=-2ZJN¯μB2g2.To calculate the total exchange energy,Eex, stored in a magnetic material, we have to integrate density of the energy (17) over the volume of the material [3] as follows:
(20)Eex=∫V‍Eexd3r=∫V‍λM2d3r+∫V‍A∑i‍(∂xiα(r))2d3r=∫V‍λM2d3r+∫V‍A[(∂xMM)2+(∂yMM)2+(∂zMM)2]d3r.The SW can be regarded as coherent precession of the magnetization vector around its equilibrium direction. Based on this observation most SW calculation are performed in linear approximation. This approximation was already used once in our paper, it is in (10). Now we will use it again to simplify the expression (20) for the total exchange energy.


### 2.1.1.  Exchange Energy in Form I

In linear approximation the magnetization vector(21)M(r)=Mx(r)x^+My(r)y^+Mz(r)z^
can be separated into two parts: a static and dynamic components. We assume that the magnetization component along the direction of the applied magnetic field, in our case it is the z-axis, is constant in time (but can be still position dependent), and its value is close to the length of the total magnetization vector as follows:
(22)M0=M0(r)≡Mz≈|M|≈const(t).
The time-depending components of the magnetization vector; Mx and My will be denoted by mx and my, respectively. We will define the dynamic magnetization vector as m=(mx,my) being a two-dimensional vector in the plane perpendicular to the direction of the saturation magnetization.The exchange energy, (20), with the help of approximation (22) can be rewritten in the following form:
(23)Eex=∫V‍λM02d3r+∫V‍A(∇mM0)2d3r.
This consists of the formula for the exchange energy, which we will call as Form I.


### 2.1.2. Exchange Energy in Form II

In the following we will make further assumptions to obtain another expression for the exchange field. We can write that(24)(∂xmM0)2=((∂xm)M0-(∂xM0)mM02)2=(∂xm)2M02+(∂xM0)2m2-M0(∂xM0)(∂xm2)M04
were in the last component of the nominator we have used 2m·∂xm=∂x(m2). The same calculations can be applied to other components of the ∇ operator in the exchange energy functional (23). The following expression for the exchange energy can be obtained:
(25)Eex=∫V‍λM02d3r+∫V‍AM02(∇m)2d3r+∫V‍AM04(∇M0)2m2d3r-∫V‍AM03(∇M0)⋅(∇m2)d3r.
In MCs the saturation magnetization is a function of the position vector with a step increase at interfaces. For bicomponent MCs (i.e., consisting of two ferromagnetic materials: A and B), M0(r) can be defined with the help of the characteristic function S(r):
(26)M0(r)=(M0,A-M0,B)S(r)+M0,B,
were
(27)S(r)={1forrinmaterialA,0forrinmaterialB,
and M0,A andM0,B are saturation magnetizations in materials A and B, respectively. We can see that in last two terms in (25) there are derivatives of M0 with respect to the position which according to (26) are derivatives from the step function, that is, ∇M0=(M0,A-M0,B)δ(r-rinterface), where δ is the Dirac delta function, and rinterface is a position vector which define the interface. It means that these two terms are connected with the exchange energy contributed only at interfaces and which are related to the jump of the saturation magnetization value (in PWM calculations we will assume parallel magnetizations in both materials). It can be shown that these two terms result in internal magnetic field components localized on interfaces and that these components introduce singularities in the equation of motion. To avoid these singularities we neglect these two terms. (It can be shown by direct calculation of functional derivatives according to (29), that these terms introduce non-Hermitian (or non-anti-Hermitian) elements into equation of motion, that is, (m/M0)∇((2A/μ0M02)∇M0), and again its physical interpretation is questionable.)To summarize, we have derived two different formulas for the exchange energy in linear approximation, which are equivalent in the case of homogeneous material. These are(28)FormI:Eex=∫V‍λM02d3r+∫V‍A(∇mM0)2d3r,FormII:Eex=∫V‍λM02d3r+∫V‍AM02(∇m)2d3r.
In this derivation we have neglected the interface anisotropy terms [9, 11]. These effects, which can be present in real materials, have a microscopic origin and are limited to the very thin area around interfaces (one or two atomic planes). In continuous effective models such effects can be included by proper effective boundary conditions imposed on dynamical component of the magnetization vector.


## 2.1.1.  Exchange Energy in Form I

In linear approximation the magnetization vector(21)M(r)=Mx(r)x^+My(r)y^+Mz(r)z^
can be separated into two parts: a static and dynamic components. We assume that the magnetization component along the direction of the applied magnetic field, in our case it is the z-axis, is constant in time (but can be still position dependent), and its value is close to the length of the total magnetization vector as follows:
(22)M0=M0(r)≡Mz≈|M|≈const(t).
The time-depending components of the magnetization vector; Mx and My will be denoted by mx and my, respectively. We will define the dynamic magnetization vector as m=(mx,my) being a two-dimensional vector in the plane perpendicular to the direction of the saturation magnetization.The exchange energy, (20), with the help of approximation (22) can be rewritten in the following form:
(23)Eex=∫V‍λM02d3r+∫V‍A(∇mM0)2d3r.
This consists of the formula for the exchange energy, which we will call as Form I.


## 2.1.2. Exchange Energy in Form II

In the following we will make further assumptions to obtain another expression for the exchange field. We can write that(24)(∂xmM0)2=((∂xm)M0-(∂xM0)mM02)2=(∂xm)2M02+(∂xM0)2m2-M0(∂xM0)(∂xm2)M04
were in the last component of the nominator we have used 2m·∂xm=∂x(m2). The same calculations can be applied to other components of the ∇ operator in the exchange energy functional (23). The following expression for the exchange energy can be obtained:
(25)Eex=∫V‍λM02d3r+∫V‍AM02(∇m)2d3r+∫V‍AM04(∇M0)2m2d3r-∫V‍AM03(∇M0)⋅(∇m2)d3r.
In MCs the saturation magnetization is a function of the position vector with a step increase at interfaces. For bicomponent MCs (i.e., consisting of two ferromagnetic materials: A and B), M0(r) can be defined with the help of the characteristic function S(r):
(26)M0(r)=(M0,A-M0,B)S(r)+M0,B,
were
(27)S(r)={1forrinmaterialA,0forrinmaterialB,
and M0,A andM0,B are saturation magnetizations in materials A and B, respectively. We can see that in last two terms in (25) there are derivatives of M0 with respect to the position which according to (26) are derivatives from the step function, that is, ∇M0=(M0,A-M0,B)δ(r-rinterface), where δ is the Dirac delta function, and rinterface is a position vector which define the interface. It means that these two terms are connected with the exchange energy contributed only at interfaces and which are related to the jump of the saturation magnetization value (in PWM calculations we will assume parallel magnetizations in both materials). It can be shown that these two terms result in internal magnetic field components localized on interfaces and that these components introduce singularities in the equation of motion. To avoid these singularities we neglect these two terms. (It can be shown by direct calculation of functional derivatives according to (29), that these terms introduce non-Hermitian (or non-anti-Hermitian) elements into equation of motion, that is, (m/M0)∇((2A/μ0M02)∇M0), and again its physical interpretation is questionable.)To summarize, we have derived two different formulas for the exchange energy in linear approximation, which are equivalent in the case of homogeneous material. These are(28)FormI:Eex=∫V‍λM02d3r+∫V‍A(∇mM0)2d3r,FormII:Eex=∫V‍λM02d3r+∫V‍AM02(∇m)2d3r.
In this derivation we have neglected the interface anisotropy terms [9, 11]. These effects, which can be present in real materials, have a microscopic origin and are limited to the very thin area around interfaces (one or two atomic planes). In continuous effective models such effects can be included by proper effective boundary conditions imposed on dynamical component of the magnetization vector.


## 2.2. Exchange Field

Exchange field can be derived from the exchange energy functionals (28) as a first variational derivative with respect to the magnetization vector [4, 59] as follows:
(29)Hex(r)=-1μ0δEexδM=-1μ0[δEexδmx,δEexδmy,δEexδM0].
This equation is written in SI units. Those variational derivatives can be calculated from Euler formula [4] as follows:
(30)δEexδmx=∂η∂mx-∇(∂η∂(∇mx))=∂η∂mx-∑i‍∂∂xi(∂η∂(∂ximx)),
where Eex=∫V‍ηd3r and η=λM02(r)+A(r)(∇m(r)/M0(r))2 or η=λM02(r)+(A(r)M02(r))(∇m(r))2 for Form I and Form II of the exchange energy, respectively (see (28)). We will perform those calculations independently for each form of the exchange energy defined in (28). During calculations we will take into account the inhomogeneity in the material, that is, spacial dependence of material parameters: A(r) and M0(r).


### 2.2.1. Form I of the Exchange Field

For Form I of the exchange energy functionalη=λM02(r)+A(r)(∇(m(r)/M0(r)))2 and we will calculate functional derivative directly from (30). First we calculate the derivative with respect of m as follows:
(31)∂η∂mx=(∂A(r)∇(m(r)/M0(r)))2∂mx=A(r)M0(r)2∇(m(r)M0(r))⋅∇(M0(r)).
This term includes the derivative of saturation magnetization, which is a step function on the interface between two magnetic materials—(26). This part of the magnetic field is localized purely at interfaces similarly as was found already in (25). This term will introduce singularity into equation of motions, and it will be neglected.After this assumption the exchange magnetic field in the Form I will be obtained solely from the second term in (30), that is,
(32)μ0Hex=∇∂η∂∇m(r)=∇∂A(r)(∇(m(r)/M0(r)))2∂∇m(r)=∇[A(r)2(∇m(r)M0(r))⋅(∂∇(m(r)/M0(r))∂∇m(r))]=∇[2A(r)M0(r)]∇m(r)M0(r).


### 2.2.2. Form II of the Exchange Field

For exchange energy written in the Form II as defined in (28), η=λM02(r)+(A(r)/M02(r))(∇m(r))2. We can calculate functional derivatives according to (30) and write Hex in a compact form without any further approximations. In this case ∂η/∂mx=0 because the first term in η is independent on ∂ximx. So the functional derivative of Eex with respect to mx take the following form:
(33)δEexδmx=-∑i‍∂∂xi(∂(A/M02)(∇m)2∂(∂ximx))=-∑i‍∂∂xi[AM02∂A∂(∂ximx)]=-∂∂x[2AM02(∂xmx)]-∂∂y[2AM02(∂ymx)]-∂∂z[2AM02(∂zmx)],
where 𝔄 denotes [(∂xmx)2+(∂xmy)2+(∂ymx)2+(∂ymy)2+(∂zmx)2+(∂zmy)2]. The same procedure can be repeated for the y component of the magnetization. Finally the exchange field can be written in the following form, that is, Form II:
(34)μ0Hex(r)=∇[2A(r)M02(r)]∇m(r).


### 2.2.3. Summary of the Exchange Field Forms

We have shown the derivation of two different expressions for the exchange field in nonuniform ferromagnetic materials in linear approximation. These are(35)FormI:Hex(r)=∇lex,I(r)∇m~(r),wherem~(r)≡m(r)M0(r),lex,I=2A(r)μ0M0(r);FormII:Hex(r)=∇lex,II(r)∇m(r),wherelex,II=2A(r)μ0M02(r).We can also add to this list the exchange field in Form III, which is derived directly from the exchange energy functional (28) (independent of the form which will be used) under assumption of the homogeneous material, that is, when the space dependence of A and M0 is not taken into account during calculations of a functional derivative (30)
(36)FormIII:Hex(r)=lex,I(r)M0(r)∇2m(r).From the parameters introduced just above:lex,I and lex,II, only the second one (i.e., from Form II) has an additional physical meaning; that is, its square root defines the exchange length [3]. (In Form III the coefficient in the exchange field is the same as in Form II, but M0 was excluded from new parameter due to simplification in the latter calculus, see (42).)It is worth to note at this moment that differential operators in the definition of the exchange field, (35)-(36), work on dynamical components of the magnetization vector, m(r) in Forms II and III, while on its normalized function in Form I, m~(r). This will make different equations and should be kept in mind during the interpretation of the eigenvectors found in PWM and boundary conditions implemented in the exchange field definitions.Those three different formulas for the exchange field will be investigated for the calculation of the magnonic band structure in thin plates of 1D and 2D Mcs with PWM.


## 2.2.1. Form I of the Exchange Field

For Form I of the exchange energy functionalη=λM02(r)+A(r)(∇(m(r)/M0(r)))2 and we will calculate functional derivative directly from (30). First we calculate the derivative with respect of m as follows:
(31)∂η∂mx=(∂A(r)∇(m(r)/M0(r)))2∂mx=A(r)M0(r)2∇(m(r)M0(r))⋅∇(M0(r)).
This term includes the derivative of saturation magnetization, which is a step function on the interface between two magnetic materials—(26). This part of the magnetic field is localized purely at interfaces similarly as was found already in (25). This term will introduce singularity into equation of motions, and it will be neglected.After this assumption the exchange magnetic field in the Form I will be obtained solely from the second term in (30), that is,
(32)μ0Hex=∇∂η∂∇m(r)=∇∂A(r)(∇(m(r)/M0(r)))2∂∇m(r)=∇[A(r)2(∇m(r)M0(r))⋅(∂∇(m(r)/M0(r))∂∇m(r))]=∇[2A(r)M0(r)]∇m(r)M0(r).


## 2.2.2. Form II of the Exchange Field

For exchange energy written in the Form II as defined in (28), η=λM02(r)+(A(r)/M02(r))(∇m(r))2. We can calculate functional derivatives according to (30) and write Hex in a compact form without any further approximations. In this case ∂η/∂mx=0 because the first term in η is independent on ∂ximx. So the functional derivative of Eex with respect to mx take the following form:
(33)δEexδmx=-∑i‍∂∂xi(∂(A/M02)(∇m)2∂(∂ximx))=-∑i‍∂∂xi[AM02∂A∂(∂ximx)]=-∂∂x[2AM02(∂xmx)]-∂∂y[2AM02(∂ymx)]-∂∂z[2AM02(∂zmx)],
where 𝔄 denotes [(∂xmx)2+(∂xmy)2+(∂ymx)2+(∂ymy)2+(∂zmx)2+(∂zmy)2]. The same procedure can be repeated for the y component of the magnetization. Finally the exchange field can be written in the following form, that is, Form II:
(34)μ0Hex(r)=∇[2A(r)M02(r)]∇m(r).


## 2.2.3. Summary of the Exchange Field Forms

We have shown the derivation of two different expressions for the exchange field in nonuniform ferromagnetic materials in linear approximation. These are(35)FormI:Hex(r)=∇lex,I(r)∇m~(r),wherem~(r)≡m(r)M0(r),lex,I=2A(r)μ0M0(r);FormII:Hex(r)=∇lex,II(r)∇m(r),wherelex,II=2A(r)μ0M02(r).We can also add to this list the exchange field in Form III, which is derived directly from the exchange energy functional (28) (independent of the form which will be used) under assumption of the homogeneous material, that is, when the space dependence of A and M0 is not taken into account during calculations of a functional derivative (30)
(36)FormIII:Hex(r)=lex,I(r)M0(r)∇2m(r).From the parameters introduced just above:lex,I and lex,II, only the second one (i.e., from Form II) has an additional physical meaning; that is, its square root defines the exchange length [3]. (In Form III the coefficient in the exchange field is the same as in Form II, but M0 was excluded from new parameter due to simplification in the latter calculus, see (42).)It is worth to note at this moment that differential operators in the definition of the exchange field, (35)-(36), work on dynamical components of the magnetization vector, m(r) in Forms II and III, while on its normalized function in Form I, m~(r). This will make different equations and should be kept in mind during the interpretation of the eigenvectors found in PWM and boundary conditions implemented in the exchange field definitions.Those three different formulas for the exchange field will be investigated for the calculation of the magnonic band structure in thin plates of 1D and 2D Mcs with PWM.


## 3. Plane Wave Method

The PWM is a useful tool used for study systems with discrete translational symmetry, including electronic, photonic, phononic and magnonic crystals [24, 27, 28, 60–66]. This method can be applied to any type of lattice and various shapes of scattering centers. The method is being constantly improved, with its field of application extending to new problems also in magnonic field [67]. Recently, the PWM has been used for the calculation of the SW spectra of 1D and 2D MCs of finite thickness [68, 69], and the magnonic spectra of thin films of 2D antidot lattices (ADLs) based on a square lattice [70, 71]. The PWM gives also a possibility for calculations of the surface effects and defect states but this requires so-called supercell formulation. The PWM in the supercell formulations was recently used to study the surface and defect influence on magnonic spectra in 2D MCs [69, 72]. For completeness we will briefly outline the PWM and explain the approximations used in this method.We will consider slabs of 1D or 2D MCs (Figure2) where the dynamics of the magnetization vector M(r,t) can be described by the LL equation as follows:
(37)∂M(r,t)∂t=-γμ0[M(r,t)×Heff(r,t)]+ξMS[M×∂M∂t],
where γ is the gyromagnetic ratio, μ0 is permeability of vacuum; as in the case of free electrons, we will assume γμ0=2.21×105 m(A s)−1. t is a time, the last term on the right describes relaxation with dimensionless damping factor ξ. The damping will be neglected in this study, while an application of the PWM for calculation of the time life of SWs in 2D MCs can be found in [67]. Heff is an effective magnetic field, which in our study will consist of three components:
(38)Heff(r,t)=H0+Hms(r,t)+Hex(r,t).H0 is a homogeneous in space and directed along the z-axis bias magnetic field, Hex(r,t) is an exchange field; its proper definitions were derived in preceding section in (35). Hms(r,t) is the demagnetizing field. In the magnetostatic approximation (with retardation effects neglected), the demagnetizing field must fulfill the magnetostatic Maxwell’s equations [59] as follows:
(39)∇×Hms(r,t)=0;∇⋅(Hms(r,t)+M(r,t))=0.
We will calculate the demagnetizing field by decomposing this field into the static and dynamic components, Hms(r) and hms(r,t), respectively. We will assume that the static part will have values different from zero only in the direction of the external magnetic field: Hms(r)=Hms(r)z^. The time dependence of the dynamic component of the demagnetizing field has the same form as that of the dynamic component of the magnetization vector: hms(r,t)=hms(r)eiωt, ω being an angular frequency of the SW.Structure of a 1D MC (a) and 2D MC (b) considered in this manuscript. 2D MC is formed by cylindrical dots A arranged in a square lattice immersed into a ferromagnetic matrix B. The external magnetic fieldH0 is applied in the direction of the z-axis. Spin waves are assumed to form standing waves in the infinite (y,z) plane or along y axis in (b) and (a), respectively. The thickness of MCs, d, is much smaller than the lattice constant a and the diameter of dots, 2R or width of stripes aA.
(a)
1D MC
(b)
2D MCIn PWM calculations we shall consider a saturated magnetization in the whole magnonic crystal. This allows us to use linear approximation and a global coordinate system in which they- and z-axes define the plane of periodicity, and the x-axis is normal to the surface of a thin plate of the MC. In the case of linear spin waves the component of the magnetization vector parallel to the static magnetic field (in this study the static magnetic field is always assumed to be oriented along the z-axis) is constant in time, and its magnitude is much greater than that of the perpendicular components: |m(r,t)|≪Mz(r) (M(r,t)=Mz(r)z^+m(r,t)). Thus, the linear approximation, introduced in the derivation of the exchange field in the previous section, can be used again, by neglecting all terms with squared m(r,t) and hms(r,t) and assuming Mz≈M0. We will only search for solutions of the LL equation corresponding to monochromatic spin waves: m(r,t)=m(r)exp(iωt).Using the linear approximation, we derive the following system of equations formx and my (and m~x and m~y) from (37) for exchange field in various formulations defined in (35). For Form I of the exchange field we get
(40)iωγμ0m~x(y)(r)=-(+)[∇⋅lex,I∇]m~y(x)(r)+(-)m~y(x)(r)(H0+Hms,z)-(+)M0h~ms,y(x)(r),
for Form II:
(41)iωγμ0mx(y)(r)=-(+)M0[∇⋅lex,II∇]my(x)(r)+(-)my(x)(r)(H0+Hms,z)-(+)M0hms,y(x)(r)
and for Form III:
(42)iωγμ0mx(y)(r)=-(+)lex,I∇2my(x)(r)+(-)my(x)(r)(H0+Hms,z)-(+)M0hms,y(x)(r).
In (40) we introduce h~ms,y(x)(r), which is a value of the demagnetizing field normalized to the saturation magnetization, M0. The formula for normalized demagnetizing field is obtained under the same approximations as used in the derivation of the exchange field and will have the same expression as hms,y(r).In MCs the material parameters, namely,A and M0, are periodic functions of the in-plane position vector r=(y,z) for 2D MC (r=(y,0) for 1D MC), with a period equal to the lattice vector a (a) as follows:
(43)M0(r+a)=M0(r),A(r+a)=A(r).
Also parameters used in formulas for the exchange field, lex,I(r) and lex,II(r), fulfill the same relation. In MCs composed of two materials each of these material parameters can be expressed by two terms: M0,A, M0,B and AA, AB, representing its respective values in each constituent material. The lattice vector a in a square lattice is any superposition of two primitive vectors: a1=az^, a2=ay^ with integer coefficients, where a is the lattice constant (see Figure 2(b)).To solve the LL equation we will use Bloch’s theorem, which asserts that a solution of a differential equation with periodic coefficients can be represented as a product of a plane-wave envelope function and a periodic function. For dynamical components of the magnetization vector and its normalized values, those are(44)m(r)=∑G‍mq(G)ei(q+G)⋅r,m~(r)=∑G‍m~q(G)ei(q+G)⋅r,
respectively. For 2D MCs G=(Gy,Gz) denotes a reciprocal lattice vector of the structure considered; in the case of square lattice G=(2π/a)(ny,nz), ny and nz are integers. The Bloch wave vector q=(qy,qz) refers to those spin waves which according to Bloch’s theorem can be limited to the first Brillouin zone (1BZ). For 1D MCs, G=(Gy,0)=((2π/a)ny,0) and q=(qy,0).In the next step we perform the Fourier transformation to map the periodic functionsMS, lex,I, and lex,II in (40)-(41) onto the reciprocal space. The transformation formulas are as follow:
(45)M0(r)=∑G‍M0(G)eiG⋅r,lex,I(r)=∑G‍lex,I(G)eiG⋅r,lex,II(r)=∑G‍lex,II(G)eiG⋅r.
In the case of cylindrical dots in 2D MC and stripes in 1D MC, the Fourier components of the saturation magnetization MS(G) and the exchange parameters, lex,I(G), lex,II(G), can be calculated analytically. The formula for the saturation magnetization in 2D MC reads as follows:
(46)M0(G)={(M0,A-M0,B)πR2a2+M0,B,forG=0,(M0,A-M0,B)2πR2a2J1(GR)GR,forG≠0,
where J1 is a Bessel function of the first kind, R is a radius of a dot. G is the length of a reciprocal wave vector G. For stripes in 1D MC with lattice constant a, M0(G) has the following form:
(47)M0(G)={(M0,A-M0,B)aAa+MS,B,forG=0,(M0,A-M0,B)2aAasin⁡GaA/2GaA/2,forG≠0,
where aA is the width of the stripe of material A. The formulas for other periodic functions of the position vector, that is, lex,I(G) and lex,II(G) have the same form.We need formulas for the static and dynamic demagnetizing fields,Hms,z(r,x), hms,x(r,x), and hms,y(r,x), to finalize the procedure, in which an eigenvalue problem in the reciprocal space is derived from LL equation. According to the ideas presented in [73], for a slab of a 2D magnonic crystal with a uniform magnetization Maxwell’s equations can be solved in the magnetostatic approximation with appropriate electromagnetic BCs at both surfaces of the slab, that is, at x=-d/2 and x=d/2. Those BCs are a continuity of the tangential components of the magnetic field vector and a normal component of the magnetic induction vector. For the considered structure, infinite in the (y,z) plane, analytical solutions in the form of Fourier series can be obtained for both the static and dynamic demagnetizing fields [68, 70] as follows:
(48)Hms,z(r,x)=-∑G‍Ms(G)G2Gz2(1-cosh⁡⁡(|G|x)e-|G|d/2)eiG⋅r;(49)hms,y(r,x)=-∑G‍my(G)|q+G|2(qy+Gy)2×(1-cosh⁡(|q+G|x)e-|q+G|d/2)ei(q+G)⋅r;(50)hms,x(r,x)=-∑G‍mx(G)cosh⁡(|q+G|x)e-|q+G|d/2ei(q+G)⋅r.
Similarly to the static component of the demagnetizing field, we take into account dynamical magnetostatic field components which depend on the same component of the magnetization vector, only. Represented in the reciprocal space for the in-plane components, these formulas for the demagnetizing fields are x dependent, that is, vary with position across the thickness of the slab. However, when the slab is thin enough (which is the case of the discussed MC, with d=30nm), the nonuniformity of the demagnetizing fields across its thickness can be neglected, and the respective field values calculated from (48)–(50) for x=0 can be used in the PWM calculations.The substitution of the (44)–(50) into (41)-(42) leads to an algebraic eigenvalue problem with eigenvalues iω/γμ0H0 as follows:
(51)M^Λmq,Λ=iωγμ0H0mq,Λ,
where Λ takes the values I, II, or III in dependence on which form of the exchange field is derived for. The eigenvector mq,ΛT=[mx,q(G1),…,mx,q(GN),my,q(G1),…,my,q(GN)] for Λ equal II or III, and mq,IT=[m~x,q(G1),…,m~x,q(GN),m~y,q(G1),…,m~y,q(GN)] when a finite number N of reciprocal lattice vectors is used in the Fourier series (44) and (45). For each form of the exchange field, the elements of the matrix M^Λ of the eigenvalue problem (51) can be written in a block-matrix form as follows:
(52)M^Λ=(MΛxxMΛxyMΛyxMΛyy).The submatrices in (52) for the exchange field in different forms are defined as follows:
(53)M^ijxx=M^ijyy=0,M^ijxy=δij+Mi,jex,Λ+(qy+Gy,j)2H0|q+Gj|2×(1-C(q+Gj,x))MS(Gi-Gj)-(Gz,i-Gz,j)2H0|Gi-Gj|2MS(Gi-Gj)(1-C(Gi-Gj,x)),M^ijxy=-δij-Mi,jex,Λ-1H0C(q+Gj,x)MS(Gi-Gj)+(Gz,i-Gz,j)2H0|Gi-Gj|2MS(Gi-Gj)(1-C(Gi-Gj,x)),
where indexes of reciprocal lattice vectors i,j,l are integers which number reciprocal lattice vectors. The additional function used in above equations is defined as follows:
(54)C(q,x)=cosh⁡(|q|x)e-|q|d/2.
Matrix elements connected with the exchange field, Mi,jex,Λ, depend on definition used. These elements are as follow:
(55)Mi,jex,I=(q+Gj)⋅(q+Gi)H0lex,I2(Gi-Gj),Mi,jex,II=∑l‍(q+Gj)⋅(q+Gl)H0lex,II2(Gl-Gj)MS(Gi-Gl),Mi,jex,III=(q+Gj)2H0lex,I2(Gi-Gj),
for Forms I, II, and III, respectively.We solve the system of (51) by standard numerical procedures designed for solving complex matrix eigenvalue problems. All the eigenvalues found by these procedures must be tested for convergence though. A satisfactory convergence of numerical solutions of (51) for all the structures considered proves to be assured by the use of 625 and 161 reciprocal lattice vectors for 2D and 1D MC, respectively.All three forms of the exchange field were used in literature in calculations of the magnonic band structure in MCs, but to our best knowledge there is missed their detailed derivations, and in this paper we would like to fill this gap. In the first paper devoted to MCs by Vasseur et al. [27], the exchange field in the form similar to Form I was postulated. The only difference is that the eigenvector is related to the dynamical component of the magnetization vector m instead its normalized value m~. The same form of the Hex was used also in other papers, for example, [67, 72, 74–76] without the derivation of its form and with the dynamical components of the magnetization vector used instead its normalized value. Of course this different definition of the eigenvectors does not influence the dispersion of spin waves found in calculations (i.e., eigenvalues) but change only the interpretation of eigenvectors. In papers [62, 63] another definition of the exchange field, that is, Form II with additional surface term, was used, and a detailed analysis of the results obtained with both formulations was presented for infinite 2D and 3D MCs. The additional component introduced in the LL equation is related to the interface anisotropy as was pointed out in the appendix of [28]. In fact this interface term can be obtained from last two terms in the exchange field defined in (25) and neglected in this study. The strict Form II of the exchange field, for the first time probably, was used by Mills in [9] where he derived the corrected Hoffman boundary conditions. In our recent papers [28, 68, 69], we used Form II for the calculations of the SW spectra in thin films of 1D and 2D MCs. The Form III of the exchange field, that is, characteristic for uniform materials, was also used for calculations of the magnonic band structure in magnetic stripes coupled by dynamic dipole interactions, [77]. Below the influence of different definitions of the exchange field on the magnonic band structure in planar MCs will be investigated.


## 4. Magnonic Spectra versus Formulation ofHex

In this section we will present results of our calculations performed with PWM for three different forms of the exchange field as defined in (35)-(36). We will present the results for a 1D MC and a 2D MC separately.


### 4.1. 1D Magnonic Crystals

We chose for our study a 1D MC consisting of Co and permalloy (Py) stripes of equal width 250 nm. The thickness of the film is 20 nm and the length is assumed infinite. Our choice is motivated by recently published papers presenting experimental and theoretical results [68, 75, 76, 78, 79]. The dispersion relation of SWs in such crystal was measured by Brillouin Light Scattering spectroscopy and was calculated from the LL equation using the finite-element method and the PWM. Very weak magnetic field, μH0=0.001 T, is applied along the stripes.We assume values of material parameters (spontaneous magnetization and exchange constant for cobalt and permalloy) equal to those presented in the experimental paper [75]. It is for cobalt: M0,Co=1.15·106 A/m, ACo=2.88·10-11 J/m and for permalloy: M0,Py=0.658·106 A/m, APy=1.1·10-11 J/m. For the gyromagnetic ratio, we assume average value proposed in [75], γ=194.6GHz/T, that is, the same for cobalt and permalloy. Three types of lines in Figure 3 represent the results of our calculations with the PWM for three different formulations of the exchange field taken from (35)-(36).Figure 3
Magnonic band structures calculated for three different definitions of the exchange field (as defined in (35)-(36)) for a thin slab of the considered 1D MC. The structure of the MC is schematically shown in the inset; stripes of Co and Py with the same width are arranged periodically with the lattice constant 500 nm. Very weak external magnetic field, μ0H0=0.001 T, is directed along stripes. The magnetostatic interaction determine, the magnonic band structure—the definition of the Hex is unimportant.The dispersion relations calculated with the three different forms of the exchange field are overlapping almost perfectly (see Figure3). It means that the effect of exchange interactions between Co and Py is minor or even negligible for this structure. It is because the magnetostatic interactions dominate over the exchange interactions for large lattice constant.To observe effects connected with various definitions of the exchange field in nonuniform materials, the decrease of the role of magnetostatic interaction is necessary. This can be done by decreasing a lattice constant. In Figure4 we show the magnonic band structure for the 1D MC with 30 nm lattice constant and 15 nm width of Co and Py stripes. We change also the film thickness to 4 nm. The rest of parameters are the same as in previous calculation.Magnonic band structure and profiles of SWs calculated for three different definitions of the exchange field for the 1D MC with the lattice constant of 30 nm and thickness 4 nm (schematically shown in the inset)—for the case of dominating exchange interaction. (a) Magnonic band structure with the Brillouin zone border marked by dashed gray line. In (b) an amplitude of the dynamical component of the magnetization vector,mx at the border of the 1st Brillouin zone is shown. Its first derivative with respect to the y is shown in (c).
(a)
(b)
(c)The differences connected with various definition of the exchange field start to be visible already near the BZ border and for higher modes also at the BZ center. The most essential difference can be observed between SW dispersions calculated according to Form I and Forms II and III. For the Form I (dashed line), the magnonic gap is absent between 1st and 2nd band, while for the other two formulations ofHex the gap exists.In Figures4(b) and 4(c) the profiles of dynamical magnetization and its first derivatives with respect to y are shown for three lowest modes. The profiles are very similar for all three forms of the exchange field. The shift of the 3rd profile along y-axis for Form I of the exchange field with respect to other two can be observed. In the inset in Figure 4(c), we can see also a discontinuity of the ∂m/∂y for the form I and II. The biggest one is observed for SWs calculated with Form I of the exchange field.


### 4.2. 2D Magnonic Crystals

Let us consider a thin film of a 2D MC composed of ferromagnetic circular dots in the square lattice and immersed in other ferromagnetic material. First, we will calculate the SW spectrum for Co dots in Py matrix in the exchange interaction dominating regime. This is obtained by assuming small lattice constanta = 30 nm and thickness d=4nm. We chose the dot diameter of 14 nm. The magnonic band structure along y-axis (i.e., perpendicular to the direction of bias magnetic field), and profiles of spin-waves are shown in Figures 5(a) and 5(b), respectively. The bias magnetic field is directed along the z-axis, and it is strong enough to saturate the sample μ0H0=0.2 T.(a) Magnonic band structure calculated for three different definitions of the exchange field for 2D MC. The MC is composed of Co dots in square lattice and immersed in Py matrix. The film has a thickness of 4 nm, the lattice constant is 30 nm, and diameter of dots is 14 nm. (b) Modulus of the dynamical component of the magnetization vector calculated for the wave vector from the center of the BZ. The profiles shown are calculated for Form II and Form I of the exchange field.
(a)
(b)We observed similar dependences as for the case of 1D MC in Figure4. For the Forms II and III of the exchange field magnonic gaps (at least partial) exist, while for the Form I the bands overlap. We found also differences in the SW profiles, especially for modes with higher frequencies.


### 4.3. Boundary Conditions

We discussed so far PWM results for various definitions of the exchange field performed for MC consisting of two materials: Co and Py only. We showed also that the different expressions for the exchange field are important only for small lattice constant. In Figures6(a) and 6(b) we show the magnonic band structure for a 1D and 2D MC, respectively, formed by Fe and yttrium iron garnet (YIG). We chose Fe and YIG because those materials have very different magnetization and exchange constants, which are MS,Fe=1.752·106 A/m, AFe=2.1·10-11 J/m, MS,YIG=0.194·106 A/m, and AYIG=0.4·10-11 J/m. The structure of the MCs is the same as in previous studies: the lattice constant 30 nm and the film thickness 4 nm for 1D MC and a = 30 nm, R = 7 nm, thickness 4 nm for 2D MC. We found that the magnonic gap (between 1st and 2nd band) exist in all band structures, also this calculated with Form I. For the 1D MC there is a big difference between the magnonic bands and the gap width calculated with Form II of the exchange field and other two forms. It is interesting to note that at Γ point there is good agreement between various formulations of Hex. In 2D MC the magnonic band structure is more complicated, gaps are smaller than for 1D but exist for all expressions of the exchange field.Magnonic band structure calculated for various forms of the exchange field. MCs composed of Fe and YIG were studied with a periodicity in (a) 1D and (b) 2D. The lattice constant is 30 nm and the film thickness 4 nm in both cases. The external magnetic field is applied alongz-axis and has the value almost 0 for 1D MC and 0.1 T for 2D MC.
(a)
(b)Now we can try to answer the question for the physical reasons of different solutions found with From I, II and III of the exchange field. We have seen significantly different results obtained from Form II, and Form I for Co/Py MC, while the solutions obtained from Forms II and III are close to each other. It is different to Fe/YIG MC where the solution for Form III is much closer to the solution of Form I. Those effects shall be related to BCs for dynamic components of the magnetization vector implemented in various formulations. We can obtain the BC implemented in the differential equation, in our case LL equation, by integrating them over the interface [9, 14, 27]. Because we are focusing on the exchange BCs, it is enough for the purpose of this study to take only terms connected with the exchange field in (41)-(42). By performing those integrations and taking a limit of a zero thickness, one can obtain the conditions of continuity of m for Forms II and III, and continuity of the m~ for Form I. At this point the difference between calculations performed in [27] and with Form I appears. In [27] the continuity of m was cast on while here is continuity of its normalized value m~. Also the second BC is obtained, that is, continuity of terms proportional to the first derivative of the dynamical component of the magnetization vector with respect to the normal to the interface. One can find that this BC requires the continuity of lex,22∂nm~,lex,II2∂nm, and ∂nm for the exchange field in Forms I, II, and III, respectively. ∂n means the derivative along the direction normal to the interface. Those boundary conditions agree with the profiles and its first derivatives found for 1D MC and shown in Figures 4(b) and 4(c), respectively (see also the inset in Figure 4(c) for the 2nd mode). The continuity of the first derivative at the interface is observed only for Form III, as expected. However, to validate which form of the exchange field is proper one and which describe properly the real physical system is out of scope of this paper. In the literature various boundary conditions were used, for example, [6–14], but the discussion under possibility for implementing them in effective continuous models is only at the beginning stage, and further investigation is required.


## 4.1. 1D Magnonic Crystals

We chose for our study a 1D MC consisting of Co and permalloy (Py) stripes of equal width 250 nm. The thickness of the film is 20 nm and the length is assumed infinite. Our choice is motivated by recently published papers presenting experimental and theoretical results [68, 75, 76, 78, 79]. The dispersion relation of SWs in such crystal was measured by Brillouin Light Scattering spectroscopy and was calculated from the LL equation using the finite-element method and the PWM. Very weak magnetic field, μH0=0.001 T, is applied along the stripes.We assume values of material parameters (spontaneous magnetization and exchange constant for cobalt and permalloy) equal to those presented in the experimental paper [75]. It is for cobalt: M0,Co=1.15·106 A/m, ACo=2.88·10-11 J/m and for permalloy: M0,Py=0.658·106 A/m, APy=1.1·10-11 J/m. For the gyromagnetic ratio, we assume average value proposed in [75], γ=194.6GHz/T, that is, the same for cobalt and permalloy. Three types of lines in Figure 3 represent the results of our calculations with the PWM for three different formulations of the exchange field taken from (35)-(36).Figure 3
Magnonic band structures calculated for three different definitions of the exchange field (as defined in (35)-(36)) for a thin slab of the considered 1D MC. The structure of the MC is schematically shown in the inset; stripes of Co and Py with the same width are arranged periodically with the lattice constant 500 nm. Very weak external magnetic field, μ0H0=0.001 T, is directed along stripes. The magnetostatic interaction determine, the magnonic band structure—the definition of the Hex is unimportant.The dispersion relations calculated with the three different forms of the exchange field are overlapping almost perfectly (see Figure3). It means that the effect of exchange interactions between Co and Py is minor or even negligible for this structure. It is because the magnetostatic interactions dominate over the exchange interactions for large lattice constant.To observe effects connected with various definitions of the exchange field in nonuniform materials, the decrease of the role of magnetostatic interaction is necessary. This can be done by decreasing a lattice constant. In Figure4 we show the magnonic band structure for the 1D MC with 30 nm lattice constant and 15 nm width of Co and Py stripes. We change also the film thickness to 4 nm. The rest of parameters are the same as in previous calculation.Magnonic band structure and profiles of SWs calculated for three different definitions of the exchange field for the 1D MC with the lattice constant of 30 nm and thickness 4 nm (schematically shown in the inset)—for the case of dominating exchange interaction. (a) Magnonic band structure with the Brillouin zone border marked by dashed gray line. In (b) an amplitude of the dynamical component of the magnetization vector,mx at the border of the 1st Brillouin zone is shown. Its first derivative with respect to the y is shown in (c).
(a)
(b)
(c)The differences connected with various definition of the exchange field start to be visible already near the BZ border and for higher modes also at the BZ center. The most essential difference can be observed between SW dispersions calculated according to Form I and Forms II and III. For the Form I (dashed line), the magnonic gap is absent between 1st and 2nd band, while for the other two formulations ofHex the gap exists.In Figures4(b) and 4(c) the profiles of dynamical magnetization and its first derivatives with respect to y are shown for three lowest modes. The profiles are very similar for all three forms of the exchange field. The shift of the 3rd profile along y-axis for Form I of the exchange field with respect to other two can be observed. In the inset in Figure 4(c), we can see also a discontinuity of the ∂m/∂y for the form I and II. The biggest one is observed for SWs calculated with Form I of the exchange field.


## 4.2. 2D Magnonic Crystals

Let us consider a thin film of a 2D MC composed of ferromagnetic circular dots in the square lattice and immersed in other ferromagnetic material. First, we will calculate the SW spectrum for Co dots in Py matrix in the exchange interaction dominating regime. This is obtained by assuming small lattice constanta = 30 nm and thickness d=4nm. We chose the dot diameter of 14 nm. The magnonic band structure along y-axis (i.e., perpendicular to the direction of bias magnetic field), and profiles of spin-waves are shown in Figures 5(a) and 5(b), respectively. The bias magnetic field is directed along the z-axis, and it is strong enough to saturate the sample μ0H0=0.2 T.(a) Magnonic band structure calculated for three different definitions of the exchange field for 2D MC. The MC is composed of Co dots in square lattice and immersed in Py matrix. The film has a thickness of 4 nm, the lattice constant is 30 nm, and diameter of dots is 14 nm. (b) Modulus of the dynamical component of the magnetization vector calculated for the wave vector from the center of the BZ. The profiles shown are calculated for Form II and Form I of the exchange field.
(a)
(b)We observed similar dependences as for the case of 1D MC in Figure4. For the Forms II and III of the exchange field magnonic gaps (at least partial) exist, while for the Form I the bands overlap. We found also differences in the SW profiles, especially for modes with higher frequencies.


## 4.3. Boundary Conditions

We discussed so far PWM results for various definitions of the exchange field performed for MC consisting of two materials: Co and Py only. We showed also that the different expressions for the exchange field are important only for small lattice constant. In Figures6(a) and 6(b) we show the magnonic band structure for a 1D and 2D MC, respectively, formed by Fe and yttrium iron garnet (YIG). We chose Fe and YIG because those materials have very different magnetization and exchange constants, which are MS,Fe=1.752·106 A/m, AFe=2.1·10-11 J/m, MS,YIG=0.194·106 A/m, and AYIG=0.4·10-11 J/m. The structure of the MCs is the same as in previous studies: the lattice constant 30 nm and the film thickness 4 nm for 1D MC and a = 30 nm, R = 7 nm, thickness 4 nm for 2D MC. We found that the magnonic gap (between 1st and 2nd band) exist in all band structures, also this calculated with Form I. For the 1D MC there is a big difference between the magnonic bands and the gap width calculated with Form II of the exchange field and other two forms. It is interesting to note that at Γ point there is good agreement between various formulations of Hex. In 2D MC the magnonic band structure is more complicated, gaps are smaller than for 1D but exist for all expressions of the exchange field.Magnonic band structure calculated for various forms of the exchange field. MCs composed of Fe and YIG were studied with a periodicity in (a) 1D and (b) 2D. The lattice constant is 30 nm and the film thickness 4 nm in both cases. The external magnetic field is applied alongz-axis and has the value almost 0 for 1D MC and 0.1 T for 2D MC.
(a)
(b)Now we can try to answer the question for the physical reasons of different solutions found with From I, II and III of the exchange field. We have seen significantly different results obtained from Form II, and Form I for Co/Py MC, while the solutions obtained from Forms II and III are close to each other. It is different to Fe/YIG MC where the solution for Form III is much closer to the solution of Form I. Those effects shall be related to BCs for dynamic components of the magnetization vector implemented in various formulations. We can obtain the BC implemented in the differential equation, in our case LL equation, by integrating them over the interface [9, 14, 27]. Because we are focusing on the exchange BCs, it is enough for the purpose of this study to take only terms connected with the exchange field in (41)-(42). By performing those integrations and taking a limit of a zero thickness, one can obtain the conditions of continuity of m for Forms II and III, and continuity of the m~ for Form I. At this point the difference between calculations performed in [27] and with Form I appears. In [27] the continuity of m was cast on while here is continuity of its normalized value m~. Also the second BC is obtained, that is, continuity of terms proportional to the first derivative of the dynamical component of the magnetization vector with respect to the normal to the interface. One can find that this BC requires the continuity of lex,22∂nm~,lex,II2∂nm, and ∂nm for the exchange field in Forms I, II, and III, respectively. ∂n means the derivative along the direction normal to the interface. Those boundary conditions agree with the profiles and its first derivatives found for 1D MC and shown in Figures 4(b) and 4(c), respectively (see also the inset in Figure 4(c) for the 2nd mode). The continuity of the first derivative at the interface is observed only for Form III, as expected. However, to validate which form of the exchange field is proper one and which describe properly the real physical system is out of scope of this paper. In the literature various boundary conditions were used, for example, [6–14], but the discussion under possibility for implementing them in effective continuous models is only at the beginning stage, and further investigation is required.


## 5. Conclusions

We presented derivations of the two different expressions for the exchange field used in literature for SW calculations in magnonic crystals with pointing at the surface terms neglected in each case. We compared these formulas with the definition of the exchange field used for SW calculations in a uniform ferromagnetic material. Numerical calculations with PWM were performed to study the influence of these different expressions on the magnonic band structure and profiles of SWs in 1D and 2D planar magnonic crystals. We found that for a large lattice constant the magnonic band structure is independent of the formulation used. It is because the magnetostatic interaction dominates over the exchange one. The situation changed for small lattice constants where in dependence on the form of the exchange field used in calculations the magnonic gap can be present or absent in magnonic band structure. By numerical calculations we showed that various formulations of the exchange field have strong relation to the boundary conditions at the interfaces between two ferromagnetic materials. Further investigation is necessary to elucidate the proper form of the exchange field which fulfill the physical boundary conditions on interfaces imposed on dynamic components of the magnetization vector in magnonic crystals.



---

*Source: 764783-2012-08-07.xml*