Theory for magnetic impurity modes in two-dimensional van der Waals ferromagnetic films

A spin-wave analysis is developed to calculate the energies of the localized excitations occurring in two-dimensional ferromagnetic van der Waals monolayers when a substitutional magnetic impurity is introduced. The magnetic ions lie on a bipartite honeycomb lattice (similar to that for graphene) and the theory includes the effects of both Ising anisotropy and single-ion anisotropy to stabilize the magnetic ordering perpendicular to the atomic plane at low temperatures. A Dyson-equation formalism, together with the spin-dependent Green’s functions derived for van der Waals monolayers, is employed to evaluate the existence conditions and energies for the impurity modes, which lie above the band of spin-wave states of the pure host material. For realistic parameter values it is found that typically two impurity modes may exist, depending on the spin quantum number for the magnetic impurity atom. Numerical applications are made to CrI3 and Cr2Ge2Te6 as the host materials.


Introduction
The recent discovery that single magnetic monolayers could be extracted from several of the previously known quasi-2D magnets, which were typically known as van der Waals (vdW) magnets, has provided great impetus to the topic of twodimensional (2D) magnetism [1][2][3].This process can often be carried out simply by exfoliation or otherwise by using standard techniques such as molecular beam epitaxy and chemical vapour deposition.The long-range order in these materials has been shown, for example, by Raman scattering [2,3].There is Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.now an expanding list of these 2D vdW ferromagnets, fabricated either as single magnetic monolayers or as films with just a few layers (see, e.g.[4][5][6][7][8] for reviews).Specific examples of vdW ferromagnetic materials, which commonly have their magnetic atoms arranged on a 2D honeycomb lattice analogous to that for graphene [9], are CrI 3 , CrBr 3 , Cr 2 Ge 2 Te 6 , and Fe 3 GeTe.Typically, the magnetization is found to be in the direction perpendicular to the plane of the layers.As well as the intrinsic scientific significance of ferromagnetic vdW materials as realizations of 2D ordered magnetic systems, they are also becoming of interest for potential device applications.For example, spintronics is a fertile area, as well as devices functioning in the fields of spin valves, spin filters, high frequency memory and storage, and neuroscience, etc (see, e.g.[10][11][12][13][14][15]).
The existence of long-range magnetic ordering in 2D materials would be excluded by the Mermin-Wagner theorem [16] for an isotropic Heisenberg magnet with exchange interactions of limited range.By contrast, in vdW magnets the well-defined magnetic ordering is stabilized with a nonzero value by symmetry-lowering magnetic anisotropies, typically assumed to be Ising-type exchange anisotropy and/or uniaxial single-ion anisotropy [3,[6][7][8].
The occurrence of an ordered magnetic ground state at temperatures well below the Curie temperature T C has prompted studies of the spin-wave (SW) dynamics in vdW materials.For example, micro-Raman scattering measurements of the two SW branches in ferromagnetic CrI 3 were made by Jin et al [17].Kapoor et al [18] employed magnetic resonance techniques to study the standing SWs excited across ultrathin films of CrI 3 .The spin-dependent transport and magnetooptical behaviour in several chromium trihalides (CrI 3 , CrBr 3 , and CrCl 3 ) were reported by Kim et al [19], who included a theoretical SW analysis.Among other theoretical studies, SW calculations were developed for bilayers of CrI 3 by Zhang et al [20].A detailed study of the SWs in ferromagnetic vdW monolayers, including temperature-dependent renormalization, was given by Mkhitaryan and Ke [21].The role of antisymmetric Dzyaloshinski-Moriya exchange interactions in a stripe geometry for antiferromagnetic vdW materials was analyzed in [22].Other types of magnetic interactions in vdW ferromagnets may include the long-range dipole-dipole interactions [23].
A quite different influence on the SW spectrum may arise, in general, from any impurities present in the system.It is well known that localized SWs may occur associated with an isolated magnetic impurity, or a low concentration of random impurities, in bulk ferromagnets and antiferromagnets.The relevant Green's-function theories to describe the scattering of SWs by impurities have been described by several authors [24][25][26][27] in the context of the Heisenberg model.This has led to the prediction of 'defect modes' occurring outside the band of bulk SWs and 'resonance modes' within the SW band.Related experimental investigations have included using Raman scattering and inelastic neutron scattering (see, e.g.[28, 29] for reviews).Subsequent work on magnetic impurities included extending the theory to Heisenberg magnetic systems in surface geometries [30,31].
Our motivation for the present work is to give an analysis for the defect SW modes occurring in a vdW ferromagnetic layer when an isolated magnetic impurity is introduced into the 2D system.This involves applying a Green's-function theory in the novel honeycomb lattice to calculate the effects of SW scattering from the impurity.We focus on predicting the localized defect modes that occur at energies above the top of the SW continuum of states.It is found that the principal types of magnetic anisotropies in vdW materials (single-ion anisotropies or Ising exchange anisotropies) affect the occurrence and energies of the modes in different ways.
The outline of this paper is as follows.In section 2 we describe the spin Hamiltonian used for the 2D magnetic system, both for the pure (host) material and with a substitutional impurity being present.The spin-spin Green's functions for the pure system are evaluated, giving the energies and spectral weights of the SW excitations.Then the role of the impurity as a SW scatterer is studied in section 3 using a Dyson-equation technique.Numerical examples are presented in section 4, including applications to CrI 3 and Cr 2 Ge 2 Te 6 , and the conclusions are in section 5.

The Green's function formalism
The honeycomb geometry of the two-dimensional (2D) vdW ferromagnetic sheet, as typified by CrI 3 among other materials, is shown in figure 1(a).There are two sublattices, denoted as A and B, and the magnetic ordering at low temperatures below the Curie temperature T C is taken along the z direction, perpendicular to the plane.The reciprocal lattice is the same as that for graphene [9], and a hexagonal cell labelled with symmetry points is depicted in figure 1(b) where the 2D wave vector is k = (k x , k y ).
The Hamiltonian for a vdW ferromagnetic will be written in a standard form [3,[6][7][8]23] that includes anisotropic exchange interactions between pairs of magnetic ions (with the anisotropy arising from an Ising term) and additional uniaxial single-ion anisotropy: Here we ignore, for simplicity, other forms of anisotropy such as those arising from magnetic dipole-dipole interactions or Dzyaloshinski-Moriya interactions.
In the pure (host) material, the exchange interaction J m,n between sites labelled m and n is assumed to take the values J 1 and J 2 between nearest neighbours (distance apart a 0 ) and next-nearest neighbours (distant apart √ 3a 0 ), respectively.Here a 0 is the bond length between magnetic ions with spin quantum number S. The Ising anisotropy along the z direction is characterized by the scalar factor ρ > 1 and the single-ion anisotropy in cases where S > 1/2 has the uniaxial form with D > 0. The final term in equation (1) represents the Zeeman energy due to an applied magnetic field B 0 in the z direction, where g is the Landé factor and µ B is the Bohr magneton.The angular brackets ⟨m, n⟩ denote a summation over distinct pairs of sites m and n only.
When a magnetic impurity ion with spin denoted by S ′ is introduced substitutionally into the lattice, the exchange interactions coupling that ion to its nearest and next-nearest neighbours may have modified values denoted by J ′ 1 and J ′ 2 , respectively, with the Ising parameter being ρ ′ .Also, we allow for modified values D ′ and g ′ for the single-ion anisotropy and Landé factor.Our methodology will be first to obtain the SW energies and unperturbed Green's functions for the pure material.Then in section 3 we introduce the impurity and analyze its effects on the SWs via a Dyson matrix formalism.

Spin wave energies
We make use of the Holstein-Primakoff transformation [32] from spin operators to boson operators.For the sites i on sublattice A the definitions appropriate to low temperatures T ≪ T C are where a † i and a i are the corresponding creation and annihilation operators.There are similar definitions for the boson operators b † j and b j at sites j on sublattice B. In terms of these operators the Hamiltonian for the pure system (now written as H 0 for emphasis) becomes Ji,j where nonlinear processes involving products of more than two operators have been ignored.Using the translational invariance symmetries in the pure material, we may transform the boson operators from the site representation to a 2D wave-vector representation with k = (k x , k y ), e.g. by writing with similar transformations applied to the b j and b † j operators.For the Fourier-transformed exchange interactions, it is straightforward to show [23] that in the honeycomb lattice, where a 0 is the lattice parameter defined earlier.The Hamiltonian of the pure host material may then be expressed as where the coefficients are defined as Alternatively, the Hamiltonian may be re-written in a matrix form as which leads to the determinantal condition for the eigenvalues E, representing the SW energies, in the form Explicitly, we may solve this equation to obtain the solutions for the two branches of the SW spectrum, where The above result is consistent with calculations by [21,23].
Its general form has been confirmed experimentally by Raman scattering [17] and will be discussed in section 4. Consistent with [23], we note that the continuum of SW energy states extends from a minimum value where is a positive anisotropy-dependent term (recalling that D > 0 and ρ > 1).

Unperturbed Green's functions
Following a standard notation (see, e.g.[33][34][35]), the Green's functions that are of interest for the spin-spin correlations in the system at energy label E are ⟨⟨a k ; X⟩⟩ E and ⟨⟨b k ; X⟩⟩ E .Here the operator X successively takes the values a † k and b † k .These quantities, which will enable us to study the SW excitations of the system, satisfy the following equations of motion [34,35]: where H 0 is the Hamiltonian quoted in equation (3) for the vdW ferromagnet.
Evaluating the commutators in equations ( 12) and ( 13) and solving the four coupled equations eventually leads to expressions for the four Green's functions as As expected, these quantities have simple poles for E at the spin-wave energies.

Theory for localized impurity modes
We next generalize the above Green's function approach to include perturbatively the effects of an isolated magnetic impurity with the modified parameters specified earlier.It is useful to write the full Hamiltonian (with an impurity present) as H = H 0 + H ′ , where H 0 is the Hamiltonian for the pure system as in equation ( 3) and H ′ is the perturbation due the impurity atom.The modified parameters associated with the latter were defined in section 2. Since we are dealing with spatial localization effects of the impurity, we now choose to work in the site representation.If we take the impurity to be at a site on sublattice A labelled as 1, there are three nearest neighbours (labelled as 2-4) and six next-nearest neighbours (labelled as 5-10), giving a cluster of 10 sites that are perturbed, as indicated in figure 2. The perturbed Hamiltonian, which affects only these sites can be written compactly as where the impurity-related parameters are defined by It is clear that a cluster size of ten sites is needed in general, but a special case of interest occurs when only the nearestneighbour exchange interactions are significant.Then the final term in equation ( 16) will vanish and the cluster size reduces to four (consisting of the impurity site and its three nearest neighbours), as discussed later.
Our objective now is to find the Green's functions G m,n (E) in the site representation when the impurity is present.By analogy with the definitions adopted in the unperturbed case, these quantities are defined as ⟨⟨a m ; a † n ⟩⟩ E if the sites m and n are both on sublattice A, and as ⟨⟨a m ; b † n ⟩⟩ E if m is on sublattice A and n is on sublattice B, and so on.In a perturbative approach, it is well known that the perturbed Green's functions are related to the unperturbed Green's functions G 0 m,n (E) for the system without an impurity by This has the usual form of a Dyson equation (see [34,35]), in which V m ′ ,n ′ represents the impurity-host interaction and so is straightforwardly related to the perturbation Hamiltonian H ′ defined in equation ( 16).In a matrix notation the above result can be expressed as , which can then be formally solved as Here I is the unit matrix and V is deduced using equation ( 16) to be in the general case.In simple terms, equation (19) represents the modified SW Green's function obtained when repeated scattering events due to the interaction V are taken into account.It is evident that, in the special case of there only being nearest-neighbour exchange interactions (such that γ 2 vanishes), the matrix V reduces to being 4 × 4.
The existence, or otherwise, of the localized impurity modes can be studied by examining the energy poles of the perturbed matrix G(E).It follows from equation (19) that any new poles (other than those that are already in G 0 (E)) must come from the [I − G 0 (E)V] −1 factor and therefore are found by the condition that The matrix elements of the spatial Green's functions G 0 (E) are just the Fourier transforms, taken with respect to the 2D wave vector k, of the quantities already found in equations ( 14) and (15).Depending on which sublattices are involved for the site representation, the transformed quantities are (for example) etc, where the domain of integration is taken over the hexagonal Brillouin zone shown (see [9,36]) in figure 1 In this work we shall restrict attention to impurity modes lying with energy E above the continuum of bulk spin wave states (i.e. the defect modes), since these are more likely to be detected experimentally compared with the 'resonance' impurity modes within the band.Also, the numerical evaluation of the transformed Green's functions will be more straightforward, as discussed later, since the energy denominators will be non-vanishing at all wave vectors.The condition for this to be the case is E must be greater than the maximum SW energy E max as defined in section 2.1.
Several of the matrix elements of G (0) will be equal to one another because of the translational and rotational symmetries of the complete 2D honeycomb lattice considered here.These symmetries have consequences for the various sets of Green's functions.For example, some simple translational symmetries give rise to the relationships , and so forth.In fact, for the case of nearest-neighbour exchange for the impurity atom, the 4 × 4 matrix has only three independent elements and can be written as where we use the shorthand , and G c = G 2,3 .For the general case, there are seven independent matrix elements; and we may write where additionally we define and G g = G 5,9 .
We may now apply the determinantal condition given in equation (21).In both the nearest-neighbour (NN) and nextnearest-neighbour (NNN) cases, it is eventually found after much algebra that we may write where we define the energy-dependent factor and H(E) is a matrix that has dimension equal to either 2 or 3 in the NN or NNN cases, respectively.First we discuss the impurity modes associated with the two factors of F(E) = 0 that are linear in Green's function components.Using the definition of ϵ, this condition can be re-arranged more conveniently as For any chosen value of the energy E, the Green's functions may be evaluated numerically and the above expression yields a value of the perturbed exchange term J ′ 1 for the impurity mode to exist.
Additional impurity modes may arise from the condition det [H(E)] = 0.In the NN case, where H(E) is a 2 × 2 matrix, its matrix elements are found to be More generally, when NNN terms are included, H(E) becomes a 3 × 3 matrix with its matrix elements defined by If we set J 2 = 0 in equation ( 29) it is easily seen that the previous determinantal result is recovered.At first sight, it may seem surprising that the H(E) in our general case is given by a 3 × 3 matrix, rather than an 8 × 8 matrix.This is explained by the simple form of the scattering matrix V in equation ( 20), which has a dimension of 10 but a rank of only 5. Specifically, its form allows us to introduce a similarity transform (to a new basis) in which the last five columns of the transformed scattering matrix all have zeroes as matrix elements, thus simplifying the form of H(E).Details of this analysis are given in the appendix.

Numerical results
We now present some numerical examples of the formalism developed in the previous two sections.First, we describe some results for the SW dispersion relations, because they

3.
The SW dispersion relation for a vdW ferromagnetic layer calculated using parameters for CrI 3 (see the text).Plots are shown for the SW energy E(k) versus wave vector components kxa 0 and kya 0 from the Brillouin zone centre (Γ point) to the corresponding zone boundary point M or K. strongly influence the spatial Green's functions entering into equations ( 26)-( 29).Then, we illustrate results predicted for the impurity modes by taking parameters appropriate to the 2D vdW ferromagnetic material CrI 3 .This is an S = 3/2 magnetic compound due to the Cr 3+ ions.It has been studied by several authors, and its important parameters are known approximately (see [21] for a summary).Here we take J 1 = 2.204 meV (≃17.8 cm −1 in wavenumber units) with J 2 /J 1 = 0.162, D/J 1 = 0.091, and ρ = 1.We will ignore effects due to the much smaller third-neighbour exchange interactions.The SW energy is shown plotted versus wave vector in figure 3, where an applied magnetic field such that gµ B B 0 /J 1 = 0.05 (or roughly B 0 ∼ 0.8 T) has been included.This is done for two wave-vector directions in the Brillouin zone, where the labels for symmetry points are given in figure 1(b).Specifically, we show variations with k x from the zone centre Γ point to a zone boundary M point and variations with k y from the Γ point to a zone boundary K point.The minimum and maximum values of E(k), which both occur at the Γ point, are 0.23 J 1 and 9.23 J 1 respectively in accordance with the expressions found in section 2.1.
It is noteworthy that the anisotropy that stabilizes the 2D magnetic ordering in CrI 3 has been attributed here to the single-ion effects, while ρ = 1.This frequently-made assumption has been brought into question in recent work (see, e.g.[37]).On investigating the effects of taking ρ ̸ = 1 (while setting the single-ion anisotropy to zero), we found that this change has only a minor effect on the form of the linearized SW dispersion relations and impurity modes for CrI 3 .
With a knowledge of the SW dispersion relations for CrI 3 , we can now evaluate the Green's functions required to solve for the impurity modes using equation ( 22) and similar expressions.When the energy E is above the threshold value of E max = 9.232 J 1 for defect modes, the Green's functions can be straightforwardly obtained numerically on employing finiteelement methods to approximate the integration over the area of the hexagonal domain.For example, some values found for the seven spatial Green's functions appearing in equations ( 23) and ( 24) are shown in table 1.It is seen, in this case, that G a (E) is the largest in magnitude and that generally the values decrease in magnitude as E increases.
The choices for the substitutional impurity atom considered to be embedded in CrI 3 can be guided by the other MX 3 vdW magnetic materials that occur in a stable form.Here M is a metallic ion and X denotes Br, Cl, or I. Thorough studies of these metallic trihalides have been made by several authors including Tomar et al [38], and some choices for M (other than Cr) are V, Ni, and Mn.The M-M bond lengths are quite similar in all these cases, leading us to consider examples in this paper where the impurity atom is taken as either V (with spin S ′ = 1/2 typically), Ni (S ′ = 1), or Mn (S ′ = 5/2).
Knowing the Green's functions at any chosen value of E, we can deduce the conditions for impurity modes (if any exist) in terms of the impurity exchange J ′ 1 by equation ( 27) and det [H(E)] = 0, together with equation (28) or (29).We shall assume, for simplicity, for the impurity atom that J ′ 2 = 0, D ′ = D, and g ′ = g, but variations of these parameters can readily be taken into account if required.Then in figure 4 we show the results for the impurity modes when both NN and NNN exchange interactions are included for the host CrI 3 material.The results are presented as plots of the energy E/J 1 versus J ′ 1 /J 1 using the full theory (10 × 10 matrix theory) with NN and NNN exchange included.Within the range of parameters considered, it is found that there are two curves in figure 4 corresponding to each choice for the magnetic impurity in terms of its impurity spin quantum number S ′ (e.g.V, Ni, or Mn).First, in each case, there is a mode arising from the condition F(E) = 0, which is shown as the solid black line.It is independent of the impurity spin label S ′ , as seen from equations ( 25)- (27), and is doubly degenerate.The other mode in each case comes from the condition det [H(E)] = 0; it is nondegenerate and does depend on S ′ .Each impurity mode emerges from the horizontal line at E max ,  shown in green, and traces out an increasing J ′ 1 /J 1 value as the energy increases.This general behaviour is analogous to what is found in many other magnetic impurity problems (see, e.g.[28,29,31]).
Next we use figure 5 to highlight the effect due to the NNN exchange interactions, which were included in the previous Figure for CrI 3 .Here in figure 5 we show the effect of setting J 2 /J 1 = 0, while keeping other parameters the same as before.As in figures 4 and 5, the dimensionless energy E/J 1 is plotted versus J ′ 1 /J 1 for each material.The impurity modes occur above the bulk SW continuum line (shown in green), which scales to be approximately the same for the two materials.
We see that the modes coming from det [H(E)] = 0, represented by the dashed lines, are shifted to a relatively small extent, but the degenerate F(E) = 0 line (full curve) is perturbed to a greater extent.
For another application of the theory to a vdW ferromagnet, we next consider a monolayer Cr 2 Ge 2 Te 6 system.Its important magnetic parameters are also known approximately from previous studies (see, e.g.[21]), and here we assume S = 3/2, J 1 = 2.01 meV (≃15.6 cm −1 in wavenumber units) with J 2 /J 1 = 0.08, D/J 1 = 0.11, ρ = 1, and gµ B B 0 /J 1 = 0.05.Compared with CrI 3 , the effects of NNN exchange interactions are weaker while the anisotropy effects are slightly greater.When the SW dispersion relation (plotted as the dimensionless quantity E(k)/J 1 versus k x a 0 or k y a 0 ) is calculated for Cr 2 Ge 2 Te 6 , the results are found to be broadly similar to figure 3 for CrI 3 .The minima and maxima for the dimensionless E(k)/J 1 are slightly different from previously; they occur when k = 0 and have the values 0.27 and 9.27, respectively, in our scaled units (compared with 0.23 and 9.23 as mentioned earlier for CrI 3 ).In meV units, using the respective J 1 for each material, we have 0.46 to 18.55 (for Cr 2 Ge 2 Te 6 ) and 0.51 to 20.34 (for CrI 3 ) as the ranges for the bulk bands in the two host materials The spatial Green's functions may be recalculated for the Cr 2 Ge 2 Te 6 monolayer as before, leading to the results for the impurity modes shown in figure 6.Here both NN and NNN exchange effects are included.In this figure we are providing a direct comparison between the modes found for Cr 2 Ge 2 Te 6 (shown in red) with those for CrI 3 (in black).In terms of the reduced units (scaled with the respective J 1 values quoted earlier), the mode energies for the two materials have a qualitatively similar behaviour.Notable differences are seen, however, for the modes derived from the condition F(E) = 0 (the solid lines), which apply for all values of the impurity spin S ′ .They are significantly split in energy, mainly because the NNN interactions are relatively much weaker in Cr 2 Ge 2 Te 6 , and these modes merge with the top of the band of bulk SW excitations at different values.

Conclusions
In conclusion we have developed here a spin-wave theory for the localized impurity modes associated with a van der Waals ferromagnetic monolayer when an impurity ion is introduced substitutionally into the 2D lattice.First, the SW modes and Green's functions for the pure host material were calculated.Then a Dyson-equation formalism was employed in a perturbative approach to deduce existence conditions for the nonresonant impurity modes above the continuum of SW states and to obtain their energies.It was found in section 3 that one type of mode, corresponding to F(E) = 0, occurred with a degeneracy factor of two and was independent of the impurity spin S ′ .The other condition, det [H(E)] = 0, typically led to one additional mode with degeneracy one that depended on S ′ .These impurity modes can be regarded as precessional spinwave modes, just as for the spin waves of the pure host material.The main physical difference is that the impurity modes are spatially localized near the impurity site and its associated cluster (of 4 or 10 sites), rather than propagating through the entire crystal.
Numerical applications were made to ferromagnetic CrI 3 and Cr 2 Ge 2 Te 6 as the host materials while V (with S ′ = 1/2), Ni (S ′ = 1) and Mn (S ′ = 5/2) were suggested as potential impurity choices.Both of these van der Waals materials have been investigated previously by Raman scattering, so the relevant parameters are fairly well known.In the case of an impure material one would expect Raman scattering to be a useful investigative technique.Any impurity modes would be observed as resonances at energies above those corresponding to the pure system, by analogy with Raman studies of other types of impure magnetic materials [29].It would be interesting, if the impurity modes in van der Waals ferromagnets are studied experimentally in the future, to deduce the J ′ 1 /J 1 value from a comparison of theory with experiment.Such an analysis has been done, for example, in some 3D bulk magnetic materials, such as the rutile-structure antiferromagnet MnF 2 with Co impurities, with moderate success [28,29].It may also be of interest in future work to generalize this low-temperature theory to higher temperatures where spin-wave renormalization effects become significant.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix. Matrix similarity transform
Some further justification for the reduction of the determinant form quoted in section 3 is provided here.We start by considering equation (20) for V.In the NN case, it simplifies to a matrix of dimension 4, as mentioned earlier.It also then has a rank of 4, which means in effect that its rows and columns are linearly independent.In the NNN case, when the dimension is 10, the rank of V is only 5, from which it follows that some linear combination of row and column operations exists that can reduce the number of nonzero rows and/or columns to 5. The set of such operations can be expressed as a similarity transform in linear algebra, i.e. by defining a transformed matrix V ′ = P −1 VP where P is a nonsingular 10 × 10 matrix.

Figure 1 .
Figure 1.(a) Geometry for the 2D space lattice as represented in the xy plane of a vdW ferromagnetic film, showing the two sublattices of magnetic sites A (blue dots) and B (red dots) forming a honeycomb lattice.(b) The 2D reciprocal lattice, showing a Brillouin-zone-centre point Γ for a hexagonal cell and some equivalent high-symmetry points that are labelled as K and M. For simplicity, only the magnetic sites are shown.

Figure 2 .
Figure 2. The cluster of ten atoms in the vdW lattice when an impurity atom, shown in green and labelled as 1, is introduced on a sublattice A site.The nearest neighbours on sublattice B (red) are labelled 2-4 and the next-nearest neighbours on sublattice A (blue) are labelled 5-10.

Figure 4 .
Figure 4. Impurity modes for CrI 3 using the full theory with NN and NNN exchange interactions included.The reduced energy E/J 1 is plotted versus J ′ 1 /J 1 for different values of the impurity spin quantum number S ′ (see the text).

Figure 5 .
Figure 5.The same as in figure 4 but showing the modified impurity modes for CrI 3 when the effects of NNN exchange interactions are ignored.

Figure 6 .
Figure 6.Impurity modes for Cr 2 Ge 2 Te 6 (shown in red) compared with those for CrI 3 (shown in black).As in figures 4 and 5, the dimensionless energy E/J 1 is plotted versus J ′ 1 /J 1 for each material.The impurity modes occur above the bulk SW continuum line (shown in green), which scales to be approximately the same for the two materials.

Table 1 .
Numerical values (expressed as multiples of J −1 1 ) obtained for the Green's functions Ga(E), G b (E), Gc(E), . . .Gg(E), that are required when evaluating the impurity modes in a CrI 3 monolayer.Results for four values of the energy E above Emax are shown.