Delving into the anisotropic interlayer exchange in bilayer CrI3

Bilayer CrI3 attracted much attention due to stacking-induced switching between the layered ferromagnetic and antiferromagnetic order. This discovery brought under the spotlight the interlayer Cr–Cr exchange interaction, which despite being much weaker than the intralayer exchange, plays an important role in shaping the magnetic properties of bilayer CrI3. In this work we delve into the anisotropic part of the interlayer exchange with the aim to separate the contributions from the Dzyaloshinskii–Moriya (DMI) and the Kitaev interactions (KI). We leverage the density functional theory calculations with spin Hamiltonian modeling and develop an energy mapping procedure to assess these anisotropic interactions with μeV accuracy. After inspecting the rhombohedral and monoclinic stacking sequences of bilayer CrI3, we reveal a considerable DMI and a weak interlayer KI between the sublattices of a monoclinic bilayer. We explain the dependence of DMI and KI on the interlayer distance, stacking sequence, and the spin–orbit coupling strength, and we suggest the dominant superexchange processes at play. In addition, we demonstrate that the single-ion anisotropy in bilayer CrI3 is highly stacking-dependent, increasing by 50% from monoclinic to rhombohedral bilayer. Remarkably, our findings prove that iodines are highly efficient in mediating the DMI across the van der Waals gap, much owing to spatially extended 5p orbitals which feature strong spin–orbit coupling. Our study gives promise that the interlayer chiral control of spin textures, demonstrated in thin metallic films where the DMI is with a much longer range, can be achieved with similar efficiency in semiconducting two-dimensional van der Waals magnets.


I. INTRODUCTION
Since the long sought discovery of two-dimensional (2D) magnets finally happened with CrI 3 and Cr 2 Ge 2 Te 6 1,2 , a large stream of scientific efforts has been directed towards achieving new capabilities with magnetic van der Waals (VdW) heterostructures [3][4][5][6][7] .Given the diversity of 2D materials, there is a profusion of possible magnetic VdW heterostructures that offer endless possibilities.Yet, to find intriguing phenomena, one doesn't need to look any further from the two layers of CrI 3 .So far, the magnetic properties of bilayer CrI 3 are manipulated by electric fields 8,9 , electrostatic doping 10 , pressure 11 , and twisting 12 .In addition, theoretical studies proposed to switch the direction of magnetization in one of its layers by spin-orbit torque 13 and predicted the magnetic photogalvanic effect 14 , magnetic polarons 15 , and magnetoelectric response in bilayer CrI 3 16 .Yet, to efficiently exploit the possibilities offered by bilayer CrI 3 in new concept devices, one has to truly understand the mechanism of the exchange coupling between the layers.
Monolayer CrI 3 is composed of chromium atoms arranged in a honeycomb lattice surrounded by edge-sharing iodine octahedra.Below the Curie temperature the S = 3/2 spins on Cr atoms are parallel and the monolayer CrI 3 is a ferromagnet (FM) with an out-of-plane magnetization [17][18][19] .The magnetic anisotropy energy (MAE), that is absolutely necessary for the long-range magnetic order to persist in 2D crystals at finite temperatures 20 , emerges in CrI 3 from an interplay between the single-ion anisotropy (SIA) and the two-ion anisotropy (TIA) occurring between the nearest neighbors Cr ions 17,21 .The TIA is usually derived within the generalized Heisenberg-Kitaev model and gives rise, in addition to the conventional isotropic Heisenberg exchange, to the Kitaev exchange and to the symmetric pseudo-dipolar interaction, depending on the bond-orientation.These terms, which we will generally label as "Kitaev-like interaction" (KI), refer to the traceless symmetric part of the most general expression for bilinear spin-spin interactions 22 and cooperate with dipole-dipole interaction in shaping the total magnetic anisotropy 23 .Like SIA, the KI in CrI 3 comes mostly from the spin-orbit coupling (SOC) on iodine atoms.MAE scales with the ligand SOC strength, being instrumental to CrI 3 showing a higher Curie temperature (T C ) than the isostructural CrBr 3 and CrCl 3 , whose ligands feature much weaker SOC than iodines [24][25][26][27] .Besides SIA and KI, SOC can give rise to the Dzyaloshinskii-Moriya interaction (DMI).However, the presence of an inversion center in the nearest neighbor Cr-Cr bonds imposes that this antisymmetric part of the anisotropic exchange is exactly zero.On the other hand the DMI is allowed between the next-nearest Cr neighbors and, although being tiny, it can play an important role in gaping the magnon spectra 28 .
In addition to magnetic properties that bilayer CrI 3 inherits from its constitutive layers, the interlayer exchange has proven extremely important as it can affect the direction of layers' magnetizations.It is an order of magnitude weaker than the intralayer exchange, but the possibility to tune it via stacking alternations made it a subject of numerous studies.For example, if we adopt the stacking from bulk CrI 3 that crystallizes either in rhombohedral (the low temperature or LT phase, R 3 space group) or monoclinic lattice (the high temperature or HT phase, C/2m space group) 29 , we end up with two different bilayer structures that we refer to as the LT and the HT structure (Fig. 1a-b).Here the theoretical 10,[30][31][32][33] and experimental 34 studies agree: the LT stacking favors the FM ordering of layers' magnetizations, whereas the HT stacking leads to layered antiferromagnetic (AFM) order.However, being realized through (at least) two iodines, the interlayer Cr-Cr coupling is mostly of super-superexchange type, which makes its microscopic description complicated due to a high number of relevant hopping processes 35,36 .
In bilayer CrI 3 studies galore the interlayer exchange is most often considered isotropic.This is a reasonable assumption, given that the bilayer CrI 3 lattice is centrosymmetric and thus the DMI is forbidden, whereas the KI, if nonzero, is usually very small.However, the strict constraints imposed on DMI by the global symmetry of the structure by no means forbid the DMI to appear locally, between the specific neighbors.Moreover, if the DMI exists between the specific neighbors and the global symmetry constraints are somehow removed -the macroscopic DMI can emerge as well.New studies warm up such expectations showing that skyrmions -topologically protected particle-like spin textures that usually appear as a consequence of DMI -can be induced via moire magnetic exchange interactions in twisted bilayer CrI 3 12,37,38

A. Computational approach
We model the bilayer CrI 3 by stacking two CrI 3 layers in rhombohedral (LT) and monoclinic (HT) sequences (Fig. 1).The structural details and computational parameters of DFT calculations are given in the Section IV.In order to study the coupling between the two spins from different layers, we need to isolate them from the rest of the magnetic environment.
To solve this problem, we use the 2 × 2 supercell and replace by Al all the Cr atoms except the two inspected ones (Fig. 1c).Al, like Cr, is trivalent so it doesn't perturb much the surrounding iodine ligand field.Therefore, the two remaining Cr atoms end up embedded into nonmagnetic crystalline environment that is reminiscent of that in bilayer CrI 3 .To check the validity of the atomic replacement method, we calculated the SIA and the intralayer nearest-neighbor exchange tensor J 1 of monolayer CrI 3 and compared the results to those obtained with the reference four-state method 41 in Supplementary Information.In calculating the interlayer exchange tensors, we follow a two-step computational procedure: in the first step the electron density and the Kohn-Sham wavefunctions are obtained from noncollinear self-consistent field calculations (SCF) without SOC.The directions of spins S 1 and S 2 on Cr atoms are constrained perpendicular to one another (Fig. 1c) using the penalty functional 42 .In the second step the Kohn-Sham wavefunctions from the first step are used, the SOC is included on all the atoms and the sum of the band energies are calculated for different directions of the spin quantization axis.The spin S 2 is kept parallel to the spin quantization axis so that it rotates together with it, whereas the direction of S 1 is unaffected by this rotation (see Fig. 1c).Given that only the sum of band energies is used, our procedure can be applied to any system provided that the use of the magnetic force theorem is justified [43][44][45] .The benchmark of our method is provided in Supplementary Information, but here we want to stress that our method can save more than 50% of computational time without affecting the accuracy of the reference four-state method 41,46 .

B. Hamiltonian of two perpendicular spins
Here we derive the model describing the coupling of two perpendicular spins, that is used for mapping the differences between DFT band energies obtained for different directions of the spin quantization axis.In closer detail, we suppose that the two spins belong to different layers, but the model works perfectly fine if spins are from the same layer.
The total energy of such two-spin system is where E nm is the nonmagnetic contribution, J is an exchange tensor that couples spins S 1 and S 2 , and A i is the SIA tensor of layer i.If the layers are of the same material and magnetic atoms are Wyckoff partners sharing the same site symmetries, we have Without loss of generality let us choose the coordinate system by fixing S 1 along the z-axis and placing S 2 in the xy-plane.Now, S 1 = (0, 0, S 1 ) and S 2 = (S 2 cos(ϕ), S 2 sin(ϕ), 0) so that the spin configuration is completely determined by a single parameter -the angle ϕ between the S 2 and the x-axis.The Eq. 1 turns into where we introduced the parameters A − ≡ A yy − A xx and A + ≡ (A xy + A yx )/2 for convenience and 2 comprises all the terms that are ϕ-independent.The important terms for our discussion are the exchange interaction energy E J (ϕ) and the contribution from the single-ion anisotropy E A (ϕ).If we interchange the roles of S 1 and S 2 the Eq. 2 yields J xz and J yz .If we further change the rotation axis from z to x we obtain J xy and J yx (and J xz and J zx that we have already).Therefore, by considering the sets of spin configurations that we denote as S 1 ∥ z ↔ S 2 ⟲ z and S 1 ∥ x ↔ S 2 ⟲ x , one obtains all the off-diagonal J -matrix elements.
In general, the J tensor decomposes into the isotropic Heisenberg exchange and the anisotropic DMI (antisymmetric) and KI (symmetric) 46 , By assumption, the spins are perpendicular and the Heisenberg exchange between them vanishes.Thus, the exchange energy contains only the DMI and the KI contributions, The DMI is usually expressed as the mixed vector product, , which introduces the Dzyaloshinskii vector, For a system with a well defined symmetry, the direction of D can be (at least partially) determined with the help of Moriya rules 22 .In Supplementary Information we derive the Moriya rules in a form that is more suitable for our purposes.
Within S 1 ∥ z ↔ S 2 ⟲ z spin configurations, the total variation of the Kitaev energy . In analogy to DMI, it is convenient to introduce the vector Note that in general the matrix K, unlike D, contains the diagonal elements as well that are out of reach of this method as their calculation require the (anti)parallel spin configurations.Therefore, strictly speaking for a given reference frame the parameter K quantifies the strength of the off-diagonal part of the Kitaev interaction.
From Eq. 2, we define the SIA constant as A = ∆E A /S 2 , where ∆E A is the total variation of the SIA energy term E A (ϕ) for the rotation of spin S around a given axis.It can be easily derived that, for a fixed rotation axis, the energy E A (ϕ) reaches extrema at and its total variation is rotates around the z-axis, the in-plane SIA constant A ∥ is obtained, whereas if it rotates around the x-(or y-) axis the fitting gives the out-of-plane SIA constant, A ⊥ .To sum up all being written about this model, if one performs the DFT calculations for S 1 ∥ z ↔ S 2 ⟲ z and S 1 ∥ x ↔ S 2 ⟲ x sets of spin configurations, the D, K, A ⊥ , and A ∥ are obtained from fitting to Eq. 2. We discuss this below for the cases of two Cr atoms from different layers.

C. Interlayer coupling of perpendicular Cr magnetic moments
Bilayer CrI 3 is an illustrative example as it reveals a few different scenarios of interlayer exchange coupling and elucidates the important role of the symmetry of local environment around Cr ions.We will refer to Cr ions pair with its surrounding ligands as "the CrI n -I n Cr cluster" (n is the number of ligands that participate in the exchange path between Cr ions).We consider the nearest and the next-nearest interlayer neighbors in both LT and HT structures.The coordinate system is chosen in a way so that the Cr-Cr bond is along the x-axis, while the z-axis is perpendicular to the CrI 3 plane (see Fig. 2a).For each considered structure, we perform the symmetry analysis to double check whether the calculated off-diagonal J -matrix elements comply with the extended Moriya rules exposed in Supplementary Information.
Starting with the LT bilayer, the nearest and the next-nearest interlayer Cr neighbors are modeled by LT 1 and LT 2 structures (see Fig. 2a-b).In LT 1 the Cr-Cr bond displays inversion, threefold rotation axis parallel to the bond, and three vertical mirror planes.
Therefore, the symmetry of CrI 3 -I 3 Cr cluster in LT 1 satisfies three Moriya rules (a,c,e) and consequently the exchange matrix must be diagonal.In complete agreement our calculations give all the zeros at the off-diagonal slots of the exchange matrix (J LT1 in Eq.S1).Moving further to next-nearest neighbors, the CrI 2 -I 2 Cr cluster in LT 2 has only the spatial-inversion symmetry (Moriya rule a) and the exchange matrix is symmetric, thus forbidding the DMI.
On the other hand, no symmetry rule forbids the KI and our calculations reveal K = 32 µeV (Table I).From this example one can actually see the Kitaev interaction in action: with S 1 fixed along the z-axis and S 2 rotating around it, the energy of the system is changing according to the cardioid pattern which is the direct consequence of the KI (see the polar plot in Fig. 2b).
Moving the discussion to HT stacking, the x-axis is a twofold rotational axis for both J zx = −J xz , and J zy = J yz .Although KI is not forbidden by symmetry, our calculations reveal its total absence in HT 1 structure and a small K of 25 µeV in HT 2 .
Strikingly, both the nearest and the next-nearest interlayer neighbors in HT structure display considerable DMI in the range of 150−200 µeV (Table I).In Moriya's seminal paper, the DMI is derived from SOC and it is shown that the DM energy is linear in SOC 22 .Given that the SOC on iodine (and not on chromium) gives the major contribution to MAE of monolayer CrI 3 17,21,47 , we assume that it is also responsible for the interlayer DMI in bilayer LT 1 (Fig. 2a) 0 0 0 0 0 0 0 0 253 0 LT 2 (Fig. 2b) 0 0 0 0 32 -14 24 -16 240 1 HT 1 (Fig. 2c) 175 0 50 -168 0 0 0 0 172 19 HT 2 (Fig. 2d) 166 0 -166 -4 25 -25 0 0 173 12 HT * 2 (CrCl-ClCr) 14 0 -14 3 1 -1 0 0 163 186 CrI 3 .This assumption can be tested through ligand replacement.The SOC constant for valence electrons in solids scales as λ ∼ 1/Z 2 , where Z is the atomic number 48 .Therefore, if iodines are replaced with chlorine, by the rule of thumb the DMI would be reduced by Cl ≈ 10 times.This simple estimate works surprisingly well, as we reveal that the D is reduced 12 times when the two iodines in HT 2 cluster are replaced by chlorine (see Fig. 2d and Table I).Moreover, this example proves the local character of the interlayer DMI, showing that only the ligands in the vicinity of Cr-Cr pair play a role in mediating this anisotropic interaction.Otherwise, the DMI would not be reduced so drastically even after the I→Cl replacement in the CrI-ICr cluster because the other iodines are still present in the structure.The dependence of DMI on SOC is further discussed in Subsection II E.
To the best of our knowledge, the effect of stacking on SIA is not addressed in previous studies of bilayer or bulk CrI 3 .Most often, it is assumed that the SIA obtained from calculations on the monolayer persists in bilayer and bulk.However, SIA is a local property and is thus sensitive to changes of the spin environment, e.g. by altering the stacking sequence in bilayer.Remarkably, we obtained the increase of 50% in A ⊥ from HT to LT stacking (Table I).Given that SIA largely contributes to magnetic anisotropy, one must take this change into account when estimating the critical temperature of bilayer or bulk CrI 3 .
Further, we obtained that LT bilayer is isotropic to the in-plane spin rotations (A ∥ = 0) like in the monolayer CrI 3 27 , but not in HT bilayer where a small in-plane SIA of A ∥ ∼ 10−20 µeV is induced by stacking.This in-plane SIA is responsible for the distortion of the cardioid that corresponds to the HT 1 structure, Fig. 2c.Taking into account all being written in this Section, one should be truly careful in ascribing the monolayer magnetic properties to bilayer and bulk VdW magnets, as the interlayer interactions proved more important than was initially expected.

D. Interlayer coupling of perpendicular magnetizations
In this Subsection we move away from the two-spin systems and provide the description of coupling between the fully magnetized layers with perpendicular magnetizations.With fully magnetized layers, we assume that all the spins in the layer point to the same direction.
In general, when spins are randomly distributed, the total energy of N × N bilayer CrI 3 (that contains 2N 2 spins per layer) is a sum of the nonmagnetic energy, the contributions from the intralayer and interlayer exchange coupling of spins, and SIA contributions at each spin site, where l = 1, 2 is the layer index and i l and j l are indices numbering the spin sites.The interlayer exchange tensor J ↕ i 1 i 2 describes the coupling between the spins i 1 and i 2 from different layers, whereas the intralayer exchange tensor J ↔ i l j l describes the coupling between spins i l and j l from the same layer l.
If we assume that all the spins belonging to one layer point to the same direction, i.e. S i l ≡ S l , the Eq. 5 greatly simplifies, where E = E N /N 2 and E nm = E N,nm /N 2 are the total and the non-magnetic energy expressed per unit cell.In Eq. 6 the contributions from A and B sublattices (see Fig. 1a,b) are separated in a sense that J ↕ A 1 B 2 describes the interlayer coupling between all the spins at A 1 sites with all the spins at B 2 sites (N 2 of each).The tensor J ↔ A describes the interaction of a single spin at A site with all the other spins from the same layer.If we specify the spins S 1 = (0, 0, S 1 ) and S 2 = (S 2 cos(ϕ), S 2 sin(ϕ), 0) the Eq. 6 reads which is the same equation as Eq. 2. The parameters Q − and Q + stemming from J ↔ and A will be used solely for fitting purposes and will not be further discussed, as our focus is on J ↕ .
We calculated the off-diagonal elements of J ↕ -matrix for HT bilayer and obtained a weak interlayer KI of K = 18 µeV.On the other hand, already from the spatial-inversion symmetry of HT bilayer we know that there is no DMI as the A 1 A 2 (A 1 B 2 ) contribution to the Dzyaloshinskii vector exactly cancels that from B 1 B 2 (B 1 A 2 ).Nevertheless, this does not mean that the interlayer DMI between the sublattices is zero.To inspect this, we modeled the A 1 A 2 sublattice of the HT structure by replacing the Cr atoms of the B 1 B 2 sublattice with Al (Fig. 3a) and obtained the Dzyaloshinskii vector with a magnitude of D A 1 A 2 = 236 µeV.This is a remarkable result, as the estimated D |J| ratio is 80% (J from Ref. 31 ), much higher than the 10% threshold which is already considered promising for skyrmionics [49][50][51] .Compared with the experimental results, the D A 1 A 2 is twice higher than the interlayer DM energy reported for TbFe/Pt/Co multilayers, which is among the strongest interlayer DMI realized in experiments 40 .

E. Dependence of DMI/KI on structural transformations and SOC
As shown in Fig. 3b, the magnitude of both the DMI and KI display a fast, exponentiallike decrease when the interlayer distance L is larger than ≈ 3 Å, consistently with the expected weak interaction between the layers across the VdW gap.At shorter separation, the evolution with L suggests a more complicated situation, especially for the KI which displays a non-monotonic behavior.Assuming a superexchange-like mechanism for both interactions, the dependence on the interlayer distance can be rationalized taking into account the possible exchange paths involving different transfer integrals (hopping terms) between Cr and I ions and their expected dependence on the relative bond lengths l.One can identify three such hopping terms, corresponding to d Cr -d Cr direct hopping (t 1 ∝ l −5 ), d Cr -p I hopping (t 2 ∝ l −7/2 ) and p I -p I hopping (t 3 ∝ l −2 ), where we adopted the bond-length dependence of Harrison 52 .In perturbation theory, the exchange interaction can be expressed quite generally as a sum of terms each with the form t n F (λ, U, ∆), where F (λ, U, ∆) is a polynomial function of atomic SOC λ, charge transfer energy ∆, and on-site Coulomb interaction U and the exponent n depends on the number of hopping terms involved in the exchange path.For instance, the contribution to the exchange interaction arising from the direct d-d hopping between Cr atoms scales as t 2 1 ∝ l −10 ∼ L −10 .Similarly, a Cr-I-Cr exchange path involving a single ligand I-p intermediate state would scale as t 4 2 ∼ L −14 , while a Cr-I-I-Cr exchange path would also include a p-p transfer integral, scaling overall as t 4 2 t 2 3 ∼ L −18 .Neglecting the details of the bond geometries and the angular dependence of each transfer integral and assuming that the dependence on the interlayer distance inherits the same scaling law of hopping terms, one can derive a general functional expression for both D and K that reads where I stands for D or K and c I n are fitting parameters effectively including all the complicated dependence on interaction matrix elements.We stress the fact that different exchange paths can give rise to perturbation terms with the same scaling dependence on the interlayer distance, that are effectively regrouped in Eq. 8.As an example, the Cr-I-Cr exchange path shares the same L −14 dependence with a process where an electron is transferred from one Cr−d to the other magnetic atom through both ligands and then transferred back via a direct d-d process.
Despite the underlying crude approximations, Eq. 8 captures surprisingly well the evolution of both DMI and KI as a function of the overall interlayer distance, as shown in Fig. 3b.Remarkably, the non-monotonic dependence of KI can be quite naturally interpreted as arising from the competition of different super-superexchange processes occurring across the VdW gap.Notwithstanding the phenomenological nature of the fitting parameters, some general trends can be deduced.For both interaction terms, c I 14 is negative while c I 10 and c I 18 are positive: the largest coefficients are those involving p-p hopping terms, which dominate at short separation, while in the opposite limit the first term -showing a slower decay -prevails.
It is worth reminding that the intralayer Kitaev interaction in monolayer CrI 3 (not to be confused with K as defined in Subsection II B) has been shown to scale quadratically with ligand iodine SOC, the transition-metal SOC contribution being negligible 21 .The superexchange mechanisms leading to anisotropic exchange interactions are therefore different from those discussed in the seminal Moriya's paper 22 , where it was assumed that the spin of the magnetic atom couples to its own orbital moment.On the other hand, the microscopic mechanisms at play are analogous to those analysed for related transition-metal dihalides, sharing with trihalides similar local bonding environments and strong magnetic anisotropies 53,54 .To shed light on SOC dependence of interlayer anisotropic exchange couplings in bilayer CrI 3 , we varied the SOC constant from λ = 0 (no SOC) to λ = 2λ 0 (double the original value), as shown in the inset of Fig. 3b, and fitted the values to the power function g(λ) = g 0 λ n .The fit yielded n D = 1.03 and n K = 2.19 (see Fig. 3c), meaning that D is a linear function of λ whereas K has nearly quadratic behavior as the intralayer Kitaev interaction.The scaling analysis on both SOC and interlayer distance suggests that superexchange mechanisms are effective in CrI 3 bilayer despite the presence of a VdW gap.
To examine the influence of a stacking sequence we gradually translated the upper layer along the a 2 lattice vector, starting from the LT and ending with the HT structure (Fig. 3c).
The different behavior of D and K with stacking alteration demonstrates the importance of angles in the Cr-I-I-Cr bonds for hopping processes that are governing the interlayer interaction.Along the inspected direction of translation, K reaches it's maximum of 33 µeV on the halfway between the LT and HT, whereas the maximum of D of 236 µeV is in the HT structure.It is interesting to see that for structures that are intermediate between LT and HT the DMI and KI are the same order of magnitude.Note that in order to find the global maximum of D and K one should inspect all possible directions, which is a demanding computational task that goes out of the scope of the present work.

DFT calculations supported by Hamiltonian modeling reveal strong interlayer DMI in
sublattice of HT bilayer CrI 3 .At the microscopic scale, DMI in HT bilayer emerges between the nearest and the next-nearest interlayer Cr neighbors due to broken local spatial inversion.However, due to the global C 2 symmetry of sublattices, the contributions from A 1 A 2 and B 1 B 2 cancel each other resulting in zero net macroscopic DMI.In addition to DMI, there is an order of magnitude weaker interlayer KI that does not vanish due to symmetry.In LT structure there is no DMI as the Cr-Cr pair with their surrounding iodines are centrosymmetric.Despite the DMI in HT structure dies out at the macroscopic scale, the main result of this work is the demonstration of the ability of iodine ligands to efficiently mediate the anisotropic exchange between the magnetic layers.This ability comes from the strong SOC of spatially extended I-5p orbitals and is impaired if iodine atoms are replaced by other ligands that have weaker SOC.The developed computational procedure offers unprecedented accuracy in calculating the anisotropic exchange interactions.Being quite general, it can be applied in any magnetic material provided that the use of the magnetic force theorem is physically sound.The symmetry analysis and the benchmark with the reference four-state method firmly support the accuracy offered by our method.In addition to the detailed analysis of the anisotropic interlayer exchange, we demonstrated that SIA heavily depends on stacking.Therefore, in estimating the Curie temperature of bilayer or bulk VdW magnet, one should not use the SIA calculated for a monolayer but instead should calculate the SIA for the system of interest.
With all being said, bilayer CrI 3 is not an appropriate 2D magnetic system for the experimental demonstration of interlayer DMI.Notwithstanding, we identified all the bricks needed to build one.Instead of attempting to modify the bilayer CrI 3 , the approach that seems more promising is to build from scratch a new heterostructure by finding an appropriate 2D magnet that can complement a layer of CrI 3 .Using a different 2D magnet as second layer has a huge advantage in inducing DMI, as one doesn't need to worry about the spatial-inversion symmetry which is trivially broken by different chemical composition of the two layers.To efficiently mediate the DMI between the magnetic ions, a candidate 2D magnet should have ligands that feature strong SOC.Most importantly, contrary to CrI 3 it should show in-plane magnetic anisotropy in order to maximize the |M 1 × M 2 | product and its MAE should be strong enough to compete with the interlayer Heisenberg exchange that favors the (anti)parallel spin configuration.
Within the CrX 3 family of 2D magnets only CrCl 3 shows an easy-plane MAE 27 , but the other properties discredit it from potential candidacy.First, its MAE is extremely small and when combined with CrI 3 all the chances are that the interlayer Heisenberg exchange will prevail, directing CrCl 3 magnetization out of plane.Second, the Cr-I-Cl-Cr interaction path is far less efficient than Cr-I-I-Cr in mediating DMI, due to small SOC on Cl.We further note that due to the lattice constant mismatch between CrI 3 and CrCl 3 one would need to match 6 × 6 structure of CrI 3 with 7 × 7 structure of CrCl 3 to build a CrI 3 /CrCl 3 heterostructure, ending up with a 680 atoms in the supercell.Hence, due to high computational demands, we were not able to check these assumptions, thus leaving the search for a potential candidate for future studies.

IV. METHODS
DFT calculations are performed using the Vasp code 55 .To describe the effects of electronic exchange and correlation we used the Perdew-Burke-Ernzerhof (PBE) functional 56 .
We did not employ and effective Hubbard parameter U but we are aware that the choice of U may affect the exchange coupling, as it is the case for monolayer CrI 3 21 .The lattice constant a = 7.005 Å of monolayer CrI 3 is obtained from spin-polarized collinear DFT calculations assuming the FM ground state.The interlayer distance is set to L opt = 3.48 Å, which corresponds to the experimental interlayer distance in bulk HT structure 29 .The bilayer made by stacking two monolayers was not relaxed any further.The lattice vector along the c-axis was set to 30 Å so that the vacuum between periodic replicas along c-axis is 20 Å thick.A cutoff of 450 eV is imposed onto the plane wave basis set and the total energies are converged to the precision of 10 −9 eV/electron.The Brillouin zone of the 2 × 2 (3 × 3) supercell was sampled by Γ-centered 4 × 4 × 1 (3 × 3 × 1) k-points mesh.The results didn't change with further increase in k-point density owing to semiconducting nature of CrI 3 .The directions of magnetic moments on Cr atoms were constrained using the approach exposed in Ref. 42 .
We carefully checked whether the size of the sphere (RWIGS) used for calculating the magnetic moments on Cr atoms and the weight of the penalty functional (LAMBDA) affect the obtained J -matrix elements.In the end we used RWIGS = 1.323Å and LAMBDA = 10.
i m100 m001 C2z C2x SA,x SB,x SB,x −SA,x −SB,x SA,x SA,y SB,y −SB,y −SA,y −SB,y −SA,y SA,z SB,z −SB,z SA,z SB,z −SA,z SB,x SA,x SA,x −SB,x −SA,x SB,x SB,y SA,y −SA,y −SB,y −SA,y −SB,y SB,z SA,z −SA,z SB,z SA,z −SB,z TABLE S2: Transformation table for axial vectors SA, SB located at points A and B under point-group symmetries for the fixed point C bisecting the straight line AB.For cases c and d we choose m001 and C2z symmetry elements without loss of generality (equivalent subcases can be obtained by rotating the reference frame around the x axis).For case e, we take n = 2 for the sake of simplicity, but the result can be easily generalized.
C 2y (case d) and C nx (case e).The transformation table for the cartesian components of the axial vectors S A , S B is given as Table S2.We notice that inversion simply acts on the site indices, swapping A and B spins, while m 001 (or m 010 ) and C nx act on cartesian components withouth swapping site indices.We focus only on symmetry constraints for off-diagonal components of the tensor J AB .Indeed, each diagonal component J AB αα always complies with the considered symmetries, as these do not mix cartesian components, as clearly seen from Table S2.
a. Inversion center at C. Applying the transformation rules of Table S2, it is found that S A,α S B,β → S B,α S A,β .It follows that, to keep the interaction term invariant, J AB αβ = J AB βα .Consistently with the first Moriya's rule, the antisymmetric component (i.e., the DMI) must vanish in the presence of inversion symmetry, and the anisotropic part of the full exchange tensor must be symmetric, with no further constraints on off-diagonal components of the pseudodipolar term.Using the Dzyaloshinskii vector and the analogous vector K we introduced to denote the off-diagonal components of the symmetric anisotropic tensor, inversion implies that D = 0 but K ̸ = 0.
b. Mirror plane perpendicular to AB and passing through C. In this case one finds: implying J AB xy = −J AB yx , J AB xz = −J AB zx and J AB yz = +J AB zy .The mirror symmetry implies that two purely antisymmetric components exist, corresponding to a Dzyaloshinskii vector D perpendicular to the AB (x) axis, or equivalently lying in the mirror plane, as stated by the second Moriya's rule.Additionally, there is a purely symmetric part of the anisotropic exchange tensor coupling spin components lying in the plane perpendicular to the AB (x) axis.
c. Mirror plane including A and B. In this situation, the reflection does not swap sites and, choosing the mirror plane to be perpendicular to the axis z (and parallel to the axis x), one gets: It follows that J AB xz = −J AB xz ≡ 0 and J AB yz = −J AB yz ≡ 0, while no symmetry constraints act on the component J AB xy .The latter condition is compatible with the third Moriya's rule, stating that D is perpendicular to the mirror plane xy, but a symmetric anisotropic part is also allowed for spin components lying within the same reflection plane.
d. Two-fold rotation axis perpendicular to AB and passing through C. Transformation rules imply that, for the two-fold rotation axis parallel to the axis z (⊥ x): S A,x S B,y → +S B,x S A,y , S A,x S B,z → −S B,x S A,z , S A,y S B,z → −S B,y S A,z .
It follows that the anisotropic exchange tensor can be decomposed in a purely antisymmetric part with D perpendicular to the two-fold axis, as J AB xz = −J AB zx and J AB yz = −J AB zy (consistenly with Moriya's fourth rule), as well as a purely symmetric part for spin components lying in a plane perpendicular to the rotation axis.S3: Generalized Moriya's rules for the anisotropic part of the exchange tensor, including both symmetric and antisymmetric components, described by the K and D vectors, respectively.All symmetry elements leave the C point bisecting the AB invariant, and their subscripts ⊥ and ∥ denote their corresponding orientation with respect to the AB direction.
e. Two-fold rotation axis parallel to AB.As for the case c, any rotation around an axis parallel to the bond AB does not swap sites, and for a two-fold rotation one gets: ≡ 0, and no symmetry constraint exists on the J AB yz component.Its antisymmetric part will give rise to a Dzyaloshinskii vector parallel to the bond axis x (as stated by the fifth Moriya's rule), while the symmetric part will couple spin components lying in a plane perpendicular to it.The Moriya's rules extended also to the symmetric part of the anisotropic exchange tensor are summarized in Table S3, in terms of the Dzyaloshinskii vector D and the K vector introduced in the main text.

III. J -MATRICES OF REPORTED STRUCTURES
Here we present the J -matrices of the structures depicted in Fig. 2

FIG. 1 :
FIG. 1: Bilayer CrI 3 with LT (a) and HT (b) stacking.Atoms from the upper layer are colored by brighter nuances.The HT structure is obtained from LT when the upper layer is translated by a 2 /3 (c) One pair of Cr atoms immersed into the CrI 3 -like environment with Cr atoms replaced by Al.Spins on Cr atoms are depicted by thick green arrows.The rotation axis is parallel to S 1 and the spin quantization axis is parallel to S 2 .

CrI 2 -I 2
FIG. 2: The structures used for modeling the nearest and the next-nearest interlayer Cr neighbors in LT and HT bilayer (upper panel).The four corresponding CrI n -I n Cr clusters are depicted in the middle panel and their symmetries are denoted.Yellow discs show the polar plots of the E(ϕ)−E 0 , where E(ϕ) is the sum of band energies for the angle ϕ of the spin quantization axis, and E 0 is the minimum of band energies around the given rotation axis.In d) the ∞-shaped energy curve corresponds to the HT * 2 structure with I atoms from CrI-ICr cluster replaced by Cl.

FIG. 3 :
FIG. 3: (a) Magnitudes of the Dzyaloshinskii (D) and Kitaev (K) vectors as a function of interlayer distance for A 1 A 2 sublattice of HT structure.Dashed lines are fitted lines with parameters c D 10 = 95 eV Å−10 , c K 10 = 13 eV Å−10 , c D 14 = −5.46× 10 3 eV Å−14 , c K 10 = −0.98 × 10 3 eV Å−14 , c D 18 = 1.11 × 10 5 eV Å−18 and c K 18 = 0.19 × 10 5 eV Å−18 .In the inset the D and K dependence on SOC constant is depicted as calculated for L = 3.48 Å (dashed lines displaying the fitting power function discussed in the text).(b) The D and K dependence on the stacking sequence from LT to HT structure for L = 3.48 Å. Continuous lines are guides for the eye.

S
A,x S B,y → −S A,x S B,y , S A,x S B,z → −S A,x S B,z , S A,y S B,z → +S A,y S B,z .It follows that J AB xy = −J AB xy ≡ 0 and J AB xz = −J AB xz

For[ 1 ]
of the Main Text and of theA 1 A 2 , A 1 B 2 , B 1 A 2and B 1 B 2 sublattices of HT bilayer CrI 3 .S = 3/2 is used the energy mapping analysis.All the values are given in µeV units.HT structure the J ↕ tensor and its decomposition onto A and B sublattices T. Moriya, Phys.Rev. 120, 91 (1960), ISSN 1536-6065.
the question is whether the iodine ligands, that also have considerable SOC, can play the role of DMI-mediator in bilayer CrI 3 .Moreover, having in mind the importance of KI in monolayer CrI 3 , how strong is the interlayer KI in bilayer CrI 3 ?We present a theoretical study that combines the density functional theory (DFT) and Hamiltonian modeling in calculating the anisotropic part of the interlayer exchange in bilayer CrI 3 .The manuscript is organized as follows: we start by presenting the employed computational approach in Subsection II A and describe the model Hamiltonian that expresses the coupling between two perpendicular spins in Subsection II B. The ability of the perpendicular-spins model to capture the changes in DFT band energies is demonstrated in Subsection II C. The validity of the model is extended to describe the coupling between fully 40Speaking of interlayer DMI in general, Vedmedenko et al.39proposed the atomistic model that predicts the formation of global chiral spin textures due to interlayer DMI between the ferromagnetic layers coupled through a nonmagnetic spacer.Recently, chiral control of spin textures is experimentally achieved in ferromagnetic TbFe/Pt/Co thin films, where the out-of-plane magnetization of TbFe is DMI-coupled with the in-plane magnetization of Co40.In this multilayer system the interlayer DMI is strong because the Pt atoms carry the conductive electrons that feature strong SOC.Having this in mind, magnetized layers in Subsection II D, where we reveal a considerable DMI and an order of magnitude weaker KI between the sublattices of bilayer CrI 3 .We demonstrate how both DMI and KI depend on structural transformations and the SOC strength in Subsection II E, suggesting possible superexchange mechanisms governing these interactions. Fily, in Section III we summarize the study by proposing a 2D magnetic heterostructure that should be a suitable platform for realizing the interlayer DMI coupling of layers' magnetizations, similar to that experimentally achieved in metallic thin films.

TABLE I :
The magnitude and the components of Dzyaloshinskii and Kitaev vectors.The last two columns present the out-of-plane and the in-plane SIA constants.