Control of magnetic states and spin interactions in bilayer CrCl3 with strain and electric fields: an ab initio study

Using ab initio density functional theory, we demonstrated the possibility of controlling the magnetic ground-state properties of bilayer CrCl\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{3}$$\end{document}3 by means of mechanical strains and electric fields. In principle, we investigated the influence of these two fields on parameters describing the spin Hamiltonian of the system. The obtained results show that biaxial strains change the magnetic ground state between ferromagnetic and antiferromagnetic phases. The mechanical strain also affects the direction and amplitude of the magnetic anisotropy energy (MAE). Importantly, the direction and amplitude of the Dzyaloshinskii–Moriya vectors are also highly tunable under external strain and electric fields. The competition between nearest-neighbor exchange interactions, MAE, and Dzyaloshinskii–Moriya interactions can lead to the stabilization of various exotic spin textures and novel magnetic excitations. The high tunability of magnetic properties by external fields makes bilayer CrCl\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{3}$$\end{document}3 a promising candidate for application in the emerging field of two-dimensional quantum spintronics and magnonics.


Introduction
The discovery of long-range magnetic order at finite temperature in atomically thin layers and the high tunability of their magnetic and electronic states represents a new route toward the next generation of ultrafast and energy-efficient spintronicbased nanodevices.Two-dimensional (2D) magnetic systems provide an ideal platform for investigating exotic properties arising from the interplay among different competing electronic and magnetic phenomena [1][2][3][4][5][6][7][8][9] .The long-range magnetic order in these 2D layers is usually stabilized at finite temperature via intrinsic magnetic anisotropies and dipolar interactions [10][11][12] .Transition chromium trihalides (CrX 3 , X= Cl, Br, I) represent a family of layered magnetic van der Waals (vdW) materials with potential applications in optoelectronics 13 and quantum spintronics 14 .Bulk CrX 3 is composed of Cr atoms arranged in hexagonal lattice layers, and each hexagonal layer is sandwiched between two halide planes, where Cr 3+ cations are octahedrally coordinated by six halide X − anions.The layers are bonded via weak vdW interactions.Bulk CrI 3 and CrBr 3 are ferromagnetic (FM) below their Curie temperatures, with magnetic moments aligned perpendicular to the Cr planes, while bulk CrCl 3 is antiferromagnetic (AFM), with a small in-plane anisotropy below its Néel temperature of 14 K 15 .The most stable stacking arrangement of bulk CrX 3 at low temperature (LT) is a rhombohedral structure with R 3 space group symmetry, while at high temperature (HT), monoclinic layer stacking with C2/m space group symmetry is the stable crystalline phase 16 .
The efficient control of spin interactions and magnetic phases in 2D systems represents important issues from technological and scientific points of view.Modifications of spin interactions cannot only tune the critical temperature of these materials by tuning their magnetic anisotropies but can also lead to the modification of magnetic states 7,9 , emerging exotic magnetic textures 17,18 , and magnonic phenomena 19 on demand.
The investigation of the magnetic properties of vdW materials such as CrX 3 shows that the spin interaction strengths between intra-and interlayer interactions are crucially different 20 In fact, compared to intralayer spin interactions, interlayer spin interactions are relatively weak and become negligible beyond the nearest neighbor atomic layers.Therefore, tuning interlayer interactions should be possible by external fields.On the other hand, since the in-plane Cr-X-Cr angle (≈ 95 • ) in these materials is close to the critical Goodenough-Kanamori-Anderson angle of 90 • , the FM superexchange interaction between Cr atoms can be externally tuned and form AFM interactions 21 .Indeed, it has been shown that the magnetic properties of 2D vdW layers can be controlled by applying strain and electric fields 8,9,[22][23][24][25][26][27] .
Monolayer and bilayer of transition chromium trihalides (CrI 3 ), with a hexagonal magnetic lattice structure, were respectively among the first 2D FM and AFM systems discovered that revealed long-range magnetic orders with a critical transition temperature of about 45 K and strong out-of-plane magnetic anisotropy that opens a gap in the magnon dispersion at the Γ-point of the magnetic Brillouin zone (BZ) 8,24,28,29 .This magnonic gap is crucial for existence of long-range magnetic order at finite temperature in 2D systems.Magnon spectrum in a hexagonal lattice with FM order has a Dirac-like dispersion at K points of the magnetic BZ 30 .Recent experimental studies show that in monolayer FM CrI 3 there is a finite gap at these K points.The origin of this gap is still unclear and under intense debate 31 .A few proposals show that this gap has a topological nature arising from either Dzyaloshinskii-Moriya or Kitaev interaction [32][33][34] .Other theories relate this gap to electron correlations and hence a nontopological origin 35 .On the other hand, bilayer of CrI 3 is an A-type AFM system, in which the intralayer exchange interaction is FM but the interlayer exchange interaction is AFM 36,37 .Voltage and strain-controlled switching between AFM and FM phases as well as tuning of the electronic band gap have been reported for this structure [21][22][23][24][25][38][39][40][41][42] .
Although there are many theoretical and experimental studies on monolayer and bulk CrX 3 materials, there are few studies on their multilayer structures.Bilayers and trilayers of these materials have greater potential for applications due to their high tunability and functionality.The weak spin interactions between adjacent layers and their magnetic ground states have not been accurately explored, which results in a poor understanding of magnetic ordering in CrX 3 bilayers.There is some discrepancy between first-principles calculations and experimental measurements regarding the layer stacking and magnetic phases in bilayer CrX 3  20, 43 .Furthermore, due to the small energy difference between FM and AFM phases, in the HT phase, different ab initio calculations of magnetic states are not always compatible with each other, and the results are sensitive to the implemented methods and initial parameters 20,44,45 .
In this paper, we investigate the magnetic properties of bilayer CrCl 3 in the presence of strain and electric fields.Monolayer CrCl 3 is a candidate for 2D-XY ferromagnetism and thus the Berezinskii-Kosterlitz-Thouless (BKT) phase transition 46 .It has been shown that the magnetic anisotropy and magnetic state of monolayer CrCl 3 can be tuned by strain and electric fields 26 .To the best of our knowledge, no systematic investigations have been conducted on the control of spin interactions in bilayer CrCl 3 .Therefore, we present a comprehensive study on the magnetic properties of bilayer CrCl 3 in the presence of biaxial strain and electric fields by means of ab initio calculations and magnetic force theory (MFT).We investigate how biaxial strain and electric fields tune the intra-and interlayer exchange interactions, Dzyaloshinskii-Moriya interactions (DMIs), and magnetic anisotropy.
The paper is organized as follows.In Sec. 2, we briefly introduce our simulation methods and our minimal model for the spin Hamiltonian.The effects of strain and electric fields on the electronic and magnetic properties of bilayer CrCl 3 are presented in Sec. 3. We conclude the paper in Sec. 4.
A schematic plot of the bilayer CrCl 3 is shown in Fig. 1.In the unit cell of monolayer CrCl 3 , there are two magnetic Cr atoms, and each Cr atom is surrounded by six I atoms in accordance with the octet rule.The two adjacent layers are coupled by weak vdW forces.The basis vector a (b) is along the zigzag (armchair) direction of the honeycomb lattice of the magnetic Cr atoms (see Fig. 1).A monolayer of CrCl 3 has three atomic sublayers.In the octahedral environment, electric crystal fields split the d orbitals of the Cr atoms into a set of triply degenerate orbitals t 2g (with lower energy) and doubly degenerate orbitals e g (with higher energy).Each of the three t 2g orbitals is occupied by one electron to minimize both the orbital energy and Coulomb interaction energy and produce an atomic magnetic moment of 3 µ B .Our ab initio calculations show a magnetic moment of 3.17 µ B for Cr atoms and an induced magnetic moment of -0.05 µ B on the Cl atoms, which is in agreement with experiments 47 .The CrCl 3 bilayer can be crystallized in rhombohedral structure with an R 3 space group symmetry at LT and monoclinic structure with the C2/m space group symmetry at HT.The total energy calculation shows that the LT phase of CrCl 3 is more stable and that its energy is lower than that of the HT phase, which is in a metastable state.Therefore, in this study, we consider only the magnetic properties of the LT phase.In the LT phase, the top layer is laterally shifted by [−1/3, 1/3] in fractional coordinates with respect to the bottom layer (Fig. 1(a)).
Ab initio density functional theory (DFT) calculations were performed using the QUANTUM ESPRESSO package 48 .In all calculations, we employed the Perdew-Burke-Ernzerhof (PBE) 49 flavor for the generalized gradient exchange-correlation functional.The BZ was sampled by a 12 × 12 × 1 k-point grid mesh, and a plane-wave cutoff energy of 80 Ry was considered.To avoid any interactions between the plane images, a 25 Å vacuum was applied along the z-axis.The total ground-state energy converged to within an accuracy of 10 −10 eV.Furthermore, the lattice parameters and atomic positions were optimized until the maximum force on each atom was less than 10 −3 eV/Å.
The vdW Grimme-D2 correction 50 was used to consider the interaction between adjacent layers.Since the onsite Hubbard interaction is essential to finding the true ground state in 2D systems, we used the DFT+U method 51 .In fact, as the band structure near the Fermi level is mostly composed of the localized d orbitals of Cr atoms, we used the DFT+U method to take into account the effect of strong electron correlations.Considering the onsite Coulomb interaction of Cr-3d orbitals, U = 3 eV in all calculations.The selected value of U in our calculations is similar to an earlier choice in the previous investigation of the CrCl 3 bilayer 47 and is around the first-principles-derived values 52 .However, it should be noted that our calculation results do not depend on Hubbard U values and remain unchanged under variation in U.The robustness of the results against different U values is shown in the Supplementary Information.
Investigating the magnetic ground state, we considered three different magnetic configurations (i) FM configuration, with all magnetic moments initialized in the same direction, (ii) A-type AFM configuration, with FM order at each layer and AFM coupling between two adjacent layers, and (iii) G-type AFM state, with magnetic moments coupled antiferromagnetically at each layer and between layers (see Fig. 1 (c)).
To compute the single ion magnetic anisotropy energy (MAE), the total energy is computed by means of fully relativistic self-consistent-field DFT calculations incorporating spin-orbit coupling (SOC) and noncollinear spin-polarization effects.The single ion MAE is defined as the difference between total energies corresponding to the magnetization orientation in-plane and out-of-plane, MAE = E IN − E OUT , and computed within MFT [53][54][55] .Therefore, a negative (positive) value of MAE indicates a uniaxial hard-axis (easy-axis) magnetic anisotropy.
To extract the spin-spin interactions in the bilayer CrCl 3 system, we adopt a minimal model spin Hamiltonian for n magnetic moments arranged in hexagonal lattices within the XY plane of where S l i denotes the magnetic moment of a Cr atom in layer l and at site i, J l,l i j represents the intra-(l = l ) and inter-(l = l ) layer symmetric Heisenberg exchange coupling between different magnetic moments at sites of i and j, D ll i j denotes the intra-(l = l ) and inter-(l = l ) layer DMI vector, and K i > 0(< 0) is the single ion easy-axis (hard-axis) magnetic anisotropy energy along the z-direction.J ll i j > 0 indicates FM coupling, while J ll i j < 0 indicates AFM coupling.The direction of the DMI vector is dictated by the symmetry of the magnetic crystal, whereas its amplitude is proportional to the SOC strength 7,[56][57][58] .The strength and sign of magnetic anisotropy are also dependent on SOC.
To obtain the spin-spin interaction parameters, presented in the spin Hamiltonian model (1), we use a Green's function method with the local rigid spin rotation, treated as a perturbation 55 .In fact, the rigid spin-rotation perturbation is done within the single-particle Green's function formalism to perturb a localized spin in the DFT electronic model.Next, spin-spin interaction parameters are mapped to the electron expressions 59 .Finally, the different parameters can be decomposed as the sum of the contributions from orbital pairs in two atoms 59 .In the next section, we compute different spin-spin interaction parameters in the spin Hamiltonian model (1).

Pristine Bilayer CrCl3
First, we study the magnetic properties of a pristine CrCl 3 bilayer in the absence of any external field.To find the magnetic ground state, we compute the energy differences between the FM state and two AFM states (types A and G; see Fig. 1(c)), E AFM − E FM .From ab initio calculations, it follows that the magnetic ground state of a CrCl 3 bilayer is an FM state 47 ; see Fig.To compute the single ion MAE, the relativistic SOC has been implemented in ab initio calculations.The MAE of the CrCl 3 bilayer is determined by calculating the energy difference between the in-plane and out-of-plane magnetic configurations; MAE=E IN − E OUT .Fig. 3 shows that the MAE, and thus the magnetic anisotropy, K i , of unstrained bilayer CrCl 3 are negative, indicating that the magnetic moments of the Cr atoms are in plane 47 .
Next, we calculate the total and orbital decomposed intra-and interlayer Heisenberg exchange interactions.Fig. 4 shows that the intralayer magnetic interactions between Cr atoms are FM (consistent with the total energy analysis), with dominant values of 1.52 and 0.29 meV for the nearest neighbor (NN) and next nearest neighbor (NNN) interactions, respectively.The NN intralayer exchange interactions between Cr atoms are mediated by e g -t 2g FM interactions and e g -e g AFM interactions, while the NNN intralayer exchange interactions result mainly from e g -t 2g FM channels.The 3rd NN interaction is AFM and contributed mostly from t 2g -t 2g coupling.Due to the longer distance between Cr atoms, the intralayer interactions between neighbors further away are negligible.Within the same layer, the Cr atom has three NN, six NNN and three 3rd NN couplings.Therefore, the sum over all intralayer neighboring pairs shows that the NN exchange interaction makes the main contribution to intralayer FM coupling (Fig. 5).
The interlayer exchange interactions between two adjacent CrCl 3 layers are computed in a similar way as intralayer interactions.The LT stacking allows one NN, nine NNN, and twelve 3rd NN couplings with larger distances than the corresponding couplings in the same layer, which results in weak interlayer exchange couplings (Fig. 4).The NN and NNN exchange couplings are FM, where the dominant contribution comes from e g -t 2g couplings.Although the AFM t 2g -t 2g channels make a noticeable contribution to first neighbor interactions, they become negligible for larger distances.Therefore, our ab initio calculations show FM interlayer coupling upon LT stacking of the CrCl 3 bilayer.
Finally, we compute different components of the DMI in the unstrained CrCl 3 bilayer; see Fig.

Effect of Biaxial Strain Field on Bilayer CrCl 3
To explore the effect of strain fields on the spin interactions of bilayer CrCl 3 , we apply an in-plane biaxial strain, defined as ε = (a − a 0 )/a 0 , where a 0 and a are lattice parameters in the absence and presence of a strain field, respectively.Our calculations show that with the application of in-plane biaxial strain fields ranging from ε = −5% to ε = +5%, the interlayer distance changes from d = 6.05 to d = 5.78.Additionally, the bond angle between the magnetic and nonmagnetic ions is changed by the strain field.In the following, we show how these changes affect different spin interactions.First, we examine the magnetic ground state in the presence of various applied strain fields.Fig. 2 summarizes our ab initio calculation results.By applying a tensile strain field ε > 0, both GGA and GGA+U calculations show that the total energy of the FM configuration is lower than that of the two other AFM configurations, and with increasing tensile strain, the FM configuration becomes more stable.
In contrast, bilayer CrCl 3 shows a quantum phase transition from the FM state to the AFM state when the compressive biaxial strain field (ε < 0) is larger than a threshold.Both GGA and GGA+U calculations show that with increasing compressive strain field, the G-type AFM state is more stable than the A-type AFM state.We obtain different critical compressive strain fields for GGA and GGA+U calculations, ε ≈ −1% and ε ≈ −3%, respectively; see Fig. 2. A similar magnetic phase transition was also reported for monolayer CrCl 3 at ε = −2.5% compressive strain 21 .
To investigate the origin of this phase transition, we compute bonding lengths and bonding angles in different configurations under various strain fields.The results are summarized in the Supplementary Information.Strain-induced changes in α Cr−Cl−Cr lead to changes in t 2g -t 2g AFM and e g -t 2g FM couplings.Consequently, the phase transition from the FM to AFM state becomes possible.
We show the evolution of spin-spin interactions under the strain in Fig. 5.The values of intralayer and interlayer symmetric Heisenberg exchange coupling increase without changing the sign under tensile strains.As a result, the FM ground state remains unchanged with increasing lattice constant.Interestingly, our GGA+U calculations show that the sign of the NN intralayer symmetric Heisenberg exchange coupling is reversed above 3% compressive strain.Fig. 5 shows that the sign of the NN intralayer symmetric exchange coefficient is negative at −5%, confirming the phase transition from an FM phase to an AFM phase under compressive strain.The NN interlayer symmetric exchange becomes negligible at 5% compressive strain due to the large interlayer distance.
The effects of strain fields on the intra-and interlayer DMIs are presented in Fig. 6.The results show that the NN and NNN intralayer components of the DMI vector increase when applying a tensile strain.As shown in Fig. 6, upon 5% tensile strain, the NN intralayers D z and D x become 10 times larger than the unstrained ones.Similarly, the NN and NNN interlayer components of the DMI vector increase with increasing tensile strain strength.Our ab initio calculations show that the NNN interlayers D z and D y are more than two times larger than the unstrained interlayers.Under 5% compressive strain, the NN intralayers D x and D y become much larger than D z , while the higher-nearest-neighbor intralayers are zero.In other words, the first-nearest-neighbor intralayer DMIs increase with decreasing intra-atomic distances, whereas the interlayer DMIs go to zero.Furthermore, Fig. 3 shows that strain fields not only modify the strength of in-plane anisotropy but also change the direction of the magnetic anisotropy from in plane to out of plane.Our calculations show that by applying a tensile strain field, the in-plane magnetic anisotropy becomes negligible around ≈ +1% and becomes out of plane for larger tensile strain fields.On the other hand, the in-plane anisotropy is increased by increasing the compressive strain fields.

Effect of Electric Field on Bilayer CrCl 3
Electric field control of magnetic states and spin interactions is also an active area of research in spintronics due to its technological applications.It was shown previously that the magnetic and electric properties of the single layer of CrX 3 can be controlled by gate voltages.In this part, we investigate the effect of perpendicular electric fields on a CrCl 3 bilayer.Our results show that within our ab initio methods, the energy difference between the AFM and FM states does not change much, and thus, the FM ground state remains the ground state of the system.Fig. 7 shows that the applied electric field increases the NN and NNN interlayer exchange interactions and even makes the third NN interactions nonzero, while the intralayer exchange couplings remain nearly unchanged.These results are consistent with the robustness of the FM state based on our total energy analyses.
The effects of perpendicular electric fields on the magnitude and direction of the intralayer and interlayer DMI vectors are depicted in Fig. 8.This figure shows that both the NN inter-and intralayer components of the DMI vector change the sign and that their amplitudes increase in the presence of a perpendicular electric field.These calculations show that at a small electric field, D z is the main component, whereas the in-plane components of the DMI vector have a dominant role at E = 1.5eV .In the presence of an electric field, the magnitude of the in-plane component (D x , D y ) of the NNN intralayer DMI increases, while the out-of-plane component (D z ) decreases.The electric field has an opposite effect on the NNN interlayer DMI, and they reduce to zero with increasing electric field.Our calculations show that an electric field up to 1.52V /nm cannot change the sign of the MAE of the CrCl 3 bilayer but increases the amplitude of the in-plane magnetic anisotropy.

Summary and Concluding Remarks
Efficient control of magnetic states and different spin interactions in a magnetic system is essential to the design of multifunctional spintronic nanodevices.Recent theoretical and experimental studies on magnetic 2D vdW materials show the high flexibility of these systems.Our ab initio calculations show that bilayer CrCl 3 is an interesting system.We show that this system can have a quantum magnetic phase transition from the FM state to the AFM state upon the application of a compressive strain field.Thus, CrX3 2D systems with FM and AFM ground states form an ideal platform for investigating magnonic phenomena.
On the other hand, we show that the amplitude and direction of uniaxial magnetic anisotropy can be changed by strain fields.The thermal stability of a 2D magnetic system is directly related to a finite bandgap in its magnon spectra.By increasing the magnetic anisotropy, one can increase the magnon gap and the corresponding critical temperature.Tuning the critical temperature is important not only for designing room temperature nanodevices but also for achieving efficient and fast magnetic switching and studying the critical behavior of quantum transport properties close to the quantum and thermal phase transitions.
Finally, we find that the DMI vectors in this system can be tuned and modified by both strain and electric fields.Both the direction and amplitude of the DMI vectors are important for designing chiral magnetic states, topological magnetic textures, and topological magnonic states, as well as engineering spin-phonon interactions.We argue that using strain engineering in CrCl 3 one can observe and study exotic phenomena, such as the magnon spin Nernst effect in the AFM state 60,61 and the realization of topological edge states in the FM state 62 , in one system.Our ab initio calculations show that while symmetric exchange and uniaxial magnetic anisotropy are sensitive to strain fields and change with a perpendicular gate voltage, the DMI vectors can be tuned by both strain and electric fields.

Figure 1 .
Figure 1.(Color online) (a), (b) Arrangement of Cr atoms in the lattice of the CrCl 3 bilayer in the LT phase.The inter-and intralayer magnetic interactions between Cr atoms are indicated.The green arrow indicates the direction of lateral shift between the top and bottom layers.(c) Schematic plot of three different spin configurations of bilayer CrCl 3 .

Figure 2 .
Figure 2. (Color online) The energy difference between the FM and AFM phases for CrCl 3 as a function of strain.The results of GGA+U (GGA) calculation indicate that a magnetic phase transition from FM to G-type AFM can occur under a compressive strain larger than 3(1)% 2. The FM ground state is robust against different values of the onsite Coulomb U.For the optimized structure of the CrCl 3 bilayer, we find a lattice constant of a 0 = 5.99, a Cr-Cl bond length of L Cr−Cl = 2.35 and a bond angle of α Cr−Cl−Cr = 94.55 • (see the Supplementary Information).

Figure 3 .
Figure 3. (Color online) Change in MAE of CrCl 3 bilayer with respect to the applied strain, within GGA and GGA+U approximations.A tensile strain equal to +2% induces a phase transition to an out-of-plane easy axis in this bilayer.

6 .
The amplitudes of the different components of the DMI vectors are small due to negligible SOC coupling in CrCl 3 .Our calculations show that for both intra-and interlayer contributions, only NN and NNN DMIs can be nonzero.Contrary to the Heisenberg exchange interactions, the intra-and interlayer DMIs have similar orders of magnitude.While the x-and z-components of the NN intralayer DMI vector, D 11 x = D 22 x and D 11 z = D 22 z , are equal and finite, the y-component of the intralayer DMI vector is zero.The situation is different for the interlayer DMI, where only the z-component of the NN DMI vector is nonzero.The finite in-plane and out-of-plane DMI vectors are promising for engineering exotic textures and helicity-dependent magnonic phenomena in this bilayer.

Figure 7 .Figure 8 .
Figure 7. (Color online) Total intralayer (left) and interlayer (right) exchange parameters of bilayer CrCl 3 as a function of electric field.A positive value indicates FM contribution, and a negative value indicates AFM contribution.