Magnon dispersion in bilayers of two-dimensional ferromagnets

We determine magnon spectra of an atomic bilayer magnet with ferromagnetic intra-and both ferro-and antiferromagnetic interlayer coupling. Analytic expressions for the full magnon band of the latter case reveal that both exchange interactions govern the fundamental magnon gap. The inter-and intralayer magnetic ordering are not independent: a stronger ferromagnetic intralayer coupling effectively strengthens the antiferromagnetic interlayer coupling as we see from comparison of two bilayer systems. The trivial topology of these exchange-anisotropy spin models without spin-orbit interaction excludes a magnon thermal Hall effect.


I. INTRODUCTION
Two-dimensional van der Waals magnets (2DvdWM) [1] are a unique platform to study magnetism in 2 + ε dimensions [2][3][4].Two-dimensional order is associated with strong intrinsic thermal fluctuations [3,5] and characteristic quantum phases [3], offering a possible test bed for competing interactions, such as Heisenberg and anisotropic exchange [6] with different range, Dzyaloshinskii-Moriya interaction (DMI) [7,8] and other spin-orbit couplings [9,10], and magnetodipolar interactions [11,12] in a rich variety of elements and crystal structures.The parameters of many properties are highly tunable by electric gating [13,14] or by strain [15,16].Of particular interest is the control of the magnetic anisotropy that modulates the spin fluctuations and allows us to study crossovers between different types of spin Hamiltonians [3].2DvdWM can be stacked with themselves or other materials into multilayers [17][18][19][20] or structured into nanodevices and directly accessed by scanning probe microscopy or other surface sensitive experimental techniques [21].
In this young field, many basic questions are still open.Only recently the magnon energy dispersion has been calculated, which is essential for understanding the spin dynamics and transport [22].For compounds with a hexagonal lattice such as CrI 3 and CrBr 3 [23] as considered here, we may expect a magnon dispersion relation similar to that of the π -electron bands of graphene-a minimum at k = 0 and two degenerate Dirac points per unit cell at an intermediate energy.This was confirmed by an analytic expression for a 2DvdWM with ferromagnetic (FM) exchange interactions [24][25][26].However, bilayers with FM intra-and inter-layer Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.Funded by Bibsam.exchange interaction show characteristic differences with bilayer graphene in terms of the degeneracy and dispersion close to the Dirac points [22].To date, the magnon dispersion for bilayers with antiferromagnetic (AFM) coupling has to the best of our knowledge been computed only numerically [22,27].
Here we extend previous theories by including a more general form of the perpendicular plane magnetic anisotropy.For the bilayer with FM intra-and AFM interlayer exchange, we report analytical results for the full spectrum by a method introduced by Colpa [28].We analyze the interplay of FM intra-and AFM interlayer couplings as reflected in the fundamental gap and total energy.The analytic solutions facilitate access to nontrivial topological properties such as the magnon Hall effect.For the class of perpendicular-plane anisotropy models without magnetization texture or spin-orbit interaction the topology is trivial, however.
The manuscript is organized as follows: In Sec.II, we define the most general spin Hamiltonian of 2DvdWM.In Sec.III, we review results on magnon spectra of an FM monolayer with different types of anisotropy and a bilayer with FM intra-and interlayer coupling.In Sec.IV we present our main results, i.e., an analytic derivation of the dispersion for a bilayer with FM intra-and AFM interlayer exchange coupling.We consider first isotropic exchange coupling for different spin configurations and subsequently include perpendicular spin anisotropy.We analyze the effect of the magnetic order on the fundamental gap as well as total energy.Finally we compute the magnon Chern numbers of the energy bands.Section V summarizes our conclusions and gives an outlook.

II. THE MODEL
Our starting point is the Heisenberg Hamiltonian with anisotropic terms that for a magnetic monolayer has the form [4] Here J i j is the exchange interaction between spins that favor ferromagnetic (J i j > 0) or antiferromagnetic (J i j < 0) order of the classical ground state, respectively.Because the exchange interaction is short-ranged, that between nearest neighbors i j dominates, while more distant ones can be disregarded.A is the single-ion anisotropy perpendicular to the plane, and α parameterizes an anisotropy in the exchange interaction in a direction α.These parameters depend on the material and can be tuned externally such as by an applied magnetic field or a gate voltage.In this paper we disregard the single-ion anisotropy (A = 0) but retain the anisotropic exchange assuming out-of-plane anisotropy, z = , x = y = 0, noting that to leading order A and are equivalent.We disregard any spin-orbit interactions at this stage.

III. REVIEW OF FM MONO-AND FM BILAYERS
We first review the Holstein-Primakoff transformation, the method of choice to treat the low frequency spectrum of spin Hamiltonians, as applied to FM monolayers [6,22] with isotropic exchange interaction (Sec.III A).Afterwards, we review different types of anisotropy in the FM coupling of the monolayer [6] (Sec.III B).Finally we consider a FM bilayer for isotropic exchange coupling as well as out-of-plane anisotropy and review the dispersion (Sec.III C) [22].This section serves essentially for fixing the geometry and the notation.

A. General method
The Holstein-Primakoff (HP) transformation of the Hamiltonian (1) replaces the local spin operators S j in favor of Boson operators a j [29,30]: At low temperatures or weak excitation we may disregard all but the zeroth order in a/ √ 2s in the series expansion of the square root.A single boson excitation a + a = 1 changes the spin projection S z = h parallel to the quantization axis z and perpendicular to the plane.After subtracting the constant ground state energy, the Hamiltonian with FM exchange interaction and zero anisotropy ( = A = 0) reads Z n.n = 3 is the number of nearest neighbors of magnetic cations on a hexagonal lattice.The lattice can be spanned by a triangular Bravais lattice with a two-atomic basis (see Fig. 1).
Transformation to momentum space leads to noninteracting magnons with energies [6] Here c k = 1 + e −i k a 1 + e −i k a 2 is the structure factor of the lattice with unit cell vectors a 1 , a 2 , as depicted in Fig. 1.This dispersion is isomorphic with the π -electrons in monolayer graphene, as shown in Fig. 2 for the first Brillouin zone (BZ).It has a minimum and maximum at the point (k = 0) and two nonequivalent Dirac cones at the K and K corners at energy 6Js with conical dispersion.

B. Anisotropies
In CrI 3 [18] the magnetic anisotropy has an easy axis along ẑ, i.e., perpendicular to the plane of the material.The Hamiltonian (1) becomes is shifted by 6 s compared to the isotropic case.This shift reflects the suppression of the Goldstone mode of rotationally symmetric systems by opening a spin wave gap at the point.
In the expansion of the HP-transformation, we restricted to leading order, thereby neglecting magnon-magnon interactions that become relevant at finite temperature.A mean-field FIG. 2. Energy dispersion of an FM monolayer with isotropic exchange coupling along high symmetry directions in the first BZ.K, K are the inequivalent Dirac points.FIG. 3. Energy dispersion for an FM monolayer with easy-plane anisotropy = −J.For further explanation see the text.treatment of higher order bosonic operators renormalizes the exchange coupling constants, and thereby also the spin wave gap [6].
We model an easy-plane anisotropic FM with J > 0 and < 0 in the Hamiltonian (1).We eliminate the nonbilinear terms of the bosonic operators a k by a Bogoliubov transformation [31], which leads to quadratic forms of Bose operators assigned to at most to two sublattices with spectrum: = −J recovers the XY model with dispersion [6,25] plotted in Fig. 3.The general monolayer Hamiltonian (8) was recently studied in Ref. [26].Note that E ± is proportional to the square root of the energy in the isotropic case.The easyplane anisotropy was observed in a monolayer of CrCl 3 [2,[32][33][34], which should therefore be a good system to study phase transitions in 2D.

C. FM bilayer
For a bilayer with FM intra-and interlayer coupling J , J ⊥ > 0 and without anisotropies we arrive at the Hamiltonian where the first and second terms describe intra-and interlayer coupling, respectively.We adopt the ratio of J ⊥ = 0.26J as predicted for CrI 3 by first-principles calculations [18].We consider here AB type stacking of 2D hexagonal lattices with a lateral shift by [2/3, 1/3] unit vectors (see Fig. 4) [18], which corresponds to the FM low-temperature crystallographic phase of bulk CrI 3 [17,19,20].We chose a unit cell for a bilayer with four atoms, A atoms A1 in the bottom-layer (1) and A2 in the top-layer (2) as well as B atoms B1 and B2 (see Fig. cell.The magnon band structure consists now of four rather than two energy bands [22] E [1] which reflects the more complex unit cell.The lowest band E [1]  − is gapless at the origin because in the absence of any anisotropy the system is invariant with respect to a global spin rotation.At the Dirac points K, K , the structure factor vanishes and E [2]  + = (12sJ + 8sJ ⊥ ), E 1 = 12sJ , where E 1 is threefold degenerate.This spectrum differs from that of the π -electrons in bilayer graphene, which are twofold degenerate at the K and K points with parabolic dispersion [35].The wave functions at the Dirac points read [1] where | jα denotes the position of the site on sublattice α.
The eigenstate [1]  K ( ) ( [2]  K ( ) ) is localized to sublattice A1(B2) in layer 1 (2).The sublattices A2 and B1 are coupled by J ⊥ , which generates an in phase or acoustic mode [3] with lower energy E 1 or out-of-phase π -shifted optical mode [4] at higher energy E [2]  + . [1], [2] , [3] correspond to excitations in which the spins on the same sublattice and layer, separated by a along (1, 0) or (− 1  2 , demonstrate that these modes also solve the Landau-Lifshitz equation for coupled classical spins. A perpendicular anisotropy can be modeled by the coupling constants J (zz) , J (zz)  ⊥ and blue shifts the frequencies, E [1]  ± = 12sJ zz + 2s Moreover, the triple degeneracy at the Dirac points is reduced to a double one.

IV. BILAYER WITH FM INTRA-AND AFM INTERLAYER EXCHANGE COUPLING
We first derive results for the dispersion for a bilayer with FM intralayer and AFM interlayer isotropic exchange interactions (Sec.IV A).Subsequently, we include a perpendicularplane anisotropy and focus on the analysis of the fundamental gap at (Sec.IV B).The topology in terms of the Berry curvature is subject of Sec.IV C.

A. Isotropic exchange interaction
Several papers discuss the impact of stacking [2,[18][19][20]36] on interlayer magnetic coupling of a CrI 3 bilayer.Depending on the type of involved interlayer orbital hybridizations, the corresponding coupling of the modeling spin Hamiltonian is FM or AFM type.For AB stacking, it has been shown by density functional theory calculations [19] that both nearest neighbor (NN) and next-nearest-neighbor (NNN) interactions determine the order of the bilayer ground state: There are one NN neighbor and 16 NNN within a unit cell, the NN contributing with AFM coupling whereas the NNN contributing with FM coupling, so that in total, interlayer magnetism in AB stacking is strongly FM.As magnetic interlayer order can be tuned by application of an electrostatic gate [14] or a magnetic field [21], however, we find it instructive to discuss both types of interlayer coupling (FM and AFM) for the same type of stacking.Here we choose AB stacking for simplicity and an AFM interlayer magnetism of the bilayer, as is induced by the NN couplings of the AB stacking.
We first calculate the energy dispersion of a bilayer with isotropic exchange coupling for different spin directions and intra/interlayer coupling strengths J /J ⊥ .The Hamiltonian reads with J , J ⊥ > 0 Again, the first sum includes the three in-plane nearest neighbors of a local moment on site i, while the second sum runs over closely spaced dimers A2, B1 between the layers.When S z = s for the spins in the top layer (2), S z = −s in the bottom layer (1) minimizes the classical ground state energy E 0 .The magnons a + i , a i are the excitations.We apply the HP transformation and expand Eq. ( 19) to leading order in the magnon operators, thereby disregarding magnon-magnon interactions, which is valid at low temperatures.In a mean-field approximation, higher terms only renormalize the exchange constants [37], as confirmed by experimental work on bilayer CrI 3 [21], at a temperature T = 0.033 J. Therefore The subscripts refer to atom α ∈ {A, B} of lattice cell i in layer ν ∈ {1, 2}.The magnon Hamiltonian then reads Ĥ − E 0 = −2sJ (i,αν),( j,α ν) ,α =α (a + j,α ν a i,αν + H.c.) As common for antiferromagnetic order, the classical ground state is not an eigenstate of the Hamiltonian since We can accommodate this issue by writing the Hamiltonian in reciprocal space as [28] Ĥ where a k = (a k,A1 , a k,B1 , a k,A2 , a k,B2 ), E c a constant to be discussed later, and D is the 8 × 8 matrix in which FIG. 5. Dispersion of a bilayer with AFM interlayer and FM intra-layer coupling, isotropic exchange coupling constants and a ratio of inter-vs intralayer coupling J ⊥ = 0.26 J .We observe a gap of order J ⊥ s at the Dirac points with quadratic instead of the linear dispersion of the FM monolayer in Fig. 2.
and c k is again the structure factor of the hexagonal lattice.Kowalska's framework [31] is not applicable for four sublattices.Instead, we diagonalize the Hamiltonian by a para-unitary transformation T of operators ( a k , a + − k ) T to the bosonic operators γ k [28]: provided that D is positive definite.E c is a further constant that will be specified below.T is paraunitary in the sense that with η = diag(I 4 , −I 4 ), where I n is the unit matrix with dimension n, which ensures that the γ (+) r, k obey bosonic commutation relations.(λ 1 , ..., λ 8 ) := (ω 1 , .., ω 4 , −ω 1 , .., −ω 4 ) are the para-values of D with paravectors v i .Equation ( 31) can be written as an value problem by multiplying by η from the left Diagonalizing the non-Hermitian matrix ηD leads to a set of four positive and four negative eigenvalues ±λ i corresponding to the 2 twofold degenerate energy bands The difference in energy bands for bilayers with AFM and FM order can be traced to the matrix η.In physical terms, two AFM-coupled sublattices (A2-B1) generate two mode families that are exchanged by a π rotation of the bilayer and hence are degenerate.The additional symmetry is also responsible for the degenerate ground state of the AFM bilayer.Breaking the interlayer symmetry by perpendicular electric and magnetic fields removes the degeneracy [22].
The dispersion (33) is plotted in Fig. 5.We find a difference E between the the zero-point energy of the magnon system and the classical ground state energy see also Eq. ( 29).The first term on the right-hand side E c = E 0 /s, arises from quantum fluctuations of the z component, while the second term reflects transverse fluctuations cause.In the following we disregard these zero-point fluctuations, but recommend their study in a future project.
Around the Dirac points K, K , the dispersion can be expanded up to second order in k as The AFM coupling J ⊥ therefore opens a gap of the order sJ ⊥ at K, K , leading to a quadratic rather than the linear dispersion found for the FM monolayer, but different effective masses.This gap implies a possible nontrivial topology.However, the Chern numbers are found to be zero for each branch, which we indicate in Sec.IV C.

B. Anisotropy
Next, we introduce an out-of-plane anisotropy with J zz > J , J zz ⊥ > J ⊥ .The matrix A then reads while B is not affected.We can still derive an analytic expression for the energy dispersion 155430-5 FIG. 6. Magnon dispersion of a bilayer with AFM inter-layer and FM intralayer coupling with anisotropic exchange coupling J zz = J xx = J yy = J.Here the interlayer couplings J ⊥ = 0.26 J , J zz ⊥ = 0.56 J and intralayer coupling J zz = 1.3 J .and plot them in Fig. 6 for coupling constants J zz = 1.3J , J zz ⊥ = 0.56 J , J ⊥ = 0.26 J .Here we adopt again a ratio of 0.26 between inter-and intralayer coupling.We assume that FM and AFM ordered layers are both AB stacked and that the ratio between inter-and intra-layer coupling (0.26 for FM CrI 3 [18]) only changes sign.Actually, AFM ordered CrI 3 has both a different (AB ) stacking and the interlayer exchange is smaller with an inter/intra layer ratio of −0.018.Other constants are known for monolayer CrI 3 [6,38] and can be tuned, for example, by an electrostatic gate [14].Here we chose them to enhance the visibility of the effects in the figures.
The anisotropy blue-shifts the lower band edge ∼J s relative to the zero-point energy E 0 − Ns (12J zz + 2J zz ⊥ ) + N k=1 4 r=1 hω r and increases the gap at the Dirac points (∼J ⊥ s for the isotropic AFM bilayer) to ∼J zz ⊥ s.We now analyze the fundamental gap hω − ( k = 0) [see Eq. (36)] plotted in Fig. 7(a) as a function of the FM coupling strength J for J ⊥ = 1.0 J 0 , J zz − J = 1.0 J 0 and J zz ⊥ − J ⊥ = 0.3 J 0 .In a simple FM the gap FM ∝ s J zz − J (37) depends on J only via anisotropy.The anisotropy gap in a pure AFM, on the other hand, depends not only on the anisotropy J zz ⊥ − J ⊥ , but also on the AFM coupling strength J ⊥ [39].The increase of the intralayer FM coupling increases the gap E − ( k = 0) according to Eq. (36), which by the reduced number of thermal magnons is equivalent to an enhanced AFM coupling.
We analyze this effect by computing the gap of a hypothetical structure in which the contributions from Eq. (37) of the FM and Eq. ( 38) of the AFM coupling at k = 0 are clearly separated.The stacking of two ferromagnetic monolayers in this "bilayer (II)" is slightly shifted such that there are two AFM-coupled dimer pairs A2 − B1 and A1 − B2 with coordination number Z AF M = 0.5 [see Fig.
does not depend explicitly on J , but on J zz − J , see Fig. 7(a) (blue line) and Eq.(39).For J = 0, the gap 3J 0 s = s(J zz − J )Z FM of bilayer (I) is governed by the anisotropy of the FM intralayer exchange only, while the AFM coupling does not contribute to the gap.The gaps converge to ∼3.61J 0 s only when the FM coupling in bilayer (I) J 5J ⊥ .This result suggests that a strong FM intralayer coupling in the realistic structure (I) increases the AFM interlayer coupling, while in the limit of weak FM coupling, the AFM order of the classical GS is less stable than in bilayer (II) [see green arrows in Fig. 7(b)].This statement is corroborated by the finite-wave vector magnon dispersion E k,0 = E − ( k) − E − ( 0) as a function of the FM coupling.The zero-k-magnon is that of an interlayer AFM in its classical GS.As E k,0 measures the energy cost of exciting a finite-k-magnon, it thereby measures the AFM coupling strength.The right panel of Fig. 8 shows an E k,0 , which indeed increases with J for both bilayers (I) and (II).The left panel of Fig. 8 shows the difference E h k,0 − E r k,0 of a hypothetical and a real bilayer for different points along the − K direction in the first BZ, which decreases with increasing J , confirming that the real bilayer approaches the effective AFM coupling strength of the hypothetical bilayer for large J .This shows that in the limit of strong intralayer coupling, magnetic order no longer depends on the choice of stacking in our specific case.

C. Topology
The topology of the magnon spectrum is reflected by the Berry curvature n k = ∇ k × u n k |i∇ k |u n k of the n = ± bands (36), where u n k is the periodic (Bloch) part of the wave function [40].For a Dirac-like spectrum, the Berry curvature is large in the vicinity of the Dirac points, which dominate the topological properties [27] as illustrated by Fig. 9. Their signs are opposite at Dirac points K, K , which means that the Chern number vanishes for each band.The topology for the bilayers in the anisotropic exchange model without spinorbit interaction is therefore trivial, without protected edge states inside the gap.The thermal Hall conductivity, which is often used to probe topological properties of systems with a Dirac-like spectrum, is proportional to the product of the Bose distribution function times the Berry curvature xy;n ( k) integrated over the first BZ [41] and vanishes as well.This corresponds to the general fact that a nonvanishing thermal Hall conductivity has so far been predicted for CrI 3 monolayer systems with the anisotropy contributions to the spin Hamiltonian of the Kitaev model [26] or the DMI [25,42].More generally, Costa et al. [43] described magnons in monolayer CrI 3 by an itinerant fermion model based on first-principles calculations, thereby circumventing model assumptions for the anisotropy.They showed that the spin-orbit coupling of iodine is essential for a nontrivial topology.

V. CONCLUSIONS
We report analytical expressions for the magnon band structure of bilayers of two-dimensional ferromagnets with (anti-) ferromagnetic interlayer exchange coupling and perpendicular anisotropy, complementing previous numerical analysis [22].An analytic expression for the fundamental gap reveals AFM and FM contributions that can be modeled by an effective coordination number.As the comparison of the spectral properties between our real bilayer system and the hypothetical toy model have shown, an increasing FM coupling in the real bilayer leads effectively to a stronger AFM interlayer coupling.The spectral properties refer to the analysis of the spectral gap as well as the energy cost associated with adding an additional magnon to the system.Both results agree with respect to the effect of stronger AFM coupling.
A natural extension of the present work would be to include next-nearest-neighbor exchange interactions, which have been shown to have an impact on magnetic interlayer coupling [18] for the AB-type stacking considered in this work.We have shown that the Chern number vanishes in the exchangeanisotropy spin model considered here, so that there is no magnon thermal Hall effect in the absence of spin-orbit interaction or complex spin texture.

FIG. 1 .
FIG. 1. Direct triangular Bravais lattice with a two-atomic basis A, B (crosses).Basis vectors a 1 , a 2 span the primitive unit cell as indicated by dashed lines.Blue circles indicate lattice point and a is the lattice constant.
FIG.4.A bilayer with AB stacking as for example in bulk BiI 3 crystals.The primitive unit cell (dashed blue lines) contains four atoms, A1 (green-rimmed black dot) of bottom layer (label 1), B2 (red-green) of top layer 2 and the stacked pair of atoms B1-A2 (black-green cross) with A2 on top of B1.The basis vectors a 1 , a 2 of the bilayer-lattice are the same as for the monolayer and are shown as blue arrows.

FIG. 7 .
FIG. 7. (a) Magnon gaps for a realistic (black line) and a hypothetical bilayer (blue line) as a function of FM coupling strength J .The anisotropy is constant with J zz − J = 1.0 J 0 , J zz ⊥ − J ⊥ = 0.3 J 0 and J ⊥ = 1.0 J 0 .(b) (left) Realistic bilayer schematic with coordination numbers Z AF M = 1 and Z FM = 3.(b) (right) Hypothetical bilayer schematic with Z AF M = 0.5 and Z FM = 3.
7(b) (right)] compared to the original coordination number Z AF M = 1 for the single dimer-pair in bilayer (I) [see Fig. 7(b) (left)].The gap of this modified system

FIG. 8 .
FIG. 8. (Left) Difference E h k,0 − E r k,0 between the hypothetical and a realistic bilayer structure as a function of FM coupling as a function of (k x , 0)[ π a ] along the − K direction in the first BZ.(Right) Energy difference E − k,0 between a magnon with wave vector (1.2, 0)[ π a ] and zero wave vector in the lower band as a function of FM coupling strength J for bilayers I (green) and II (violet).

FIG. 9 .
FIG. 9.The Berry curvature xy;n ( k) for the bands E + (left) and E − (right) of a bilayer with AFM interlayer and FM intralayer coupling.The exchange coupling constants are chosen as in Fig. 6.