Moir\'e-driven multiferroic order in twisted CrCl$_3$, CrBr$_3$ and CrI$_3$ bilayers

Layered van der Waals materials have risen as a powerful platform to engineer artificial competing states of matter. Here we show the emergence of multiferroic order in twisted chromium trihalide bilayers, an order fully driven by the moir\'e pattern and absent in aligned multilayers. Using a combination of spin models and ab initio calculations, we show that a spin texture is generated in the moir\'e supercell of the twisted system as a consequence of the competition between stacking-dependent interlayer magnetic exchange and magnetic anisotropy. An electric polarization arises associated with such a non-collinear magnetic state due to the spin-orbit coupling, leading to the emergence of a local ferroelectric order following the moir\'e. Among the stochiometric trihalides, our results show that twisted CrBr$_3$ bilayers give rise to the strongest multiferroic order. We further show the emergence of a strong magnetoelectric coupling, which allows the electric generation and control of magnetic skyrmions. Our results put forward twisted chromium trihalide bilayers, and in particular CrBr$_3$ bilayers, as a powerful platform to engineer artificial multiferroic materials and electrically tunable topological magnetic textures.


I. INTRODUCTION
Multiferroic materials display more than one ferroic order at the same time 1-3 , and in particular, they can simultaneously host magnetic and ferroelectric orders. The existence of multiple symmetry-breaking orders allows having a coupling between electric and magnetic degrees of freedom 4 .
Here we demonstrate the emergence of a multiferroic state in the family of twisted chromium trihalide CrX 3 (X=Cl, Br and I) bilayers by combining first principles and effective spin Hamiltonians. We first show the emergence of a non-collinear spin texture due to the modulation of the interlayer exchange coupling in the moiré unit cell. Associated with the spin texture an electric polarization emerges as a consequence of spin-orbit coupling (SOC) and the local magnetic non-collinearity in the moiré domains. Using ab initio calculations we extract the value of the electric polarization driven by noncollinear magnetic texture, and demonstrate its dependence on the halide of CrX 3 . Therefore, we provide a quantification of the resulting multiferroic order. Furthermore, we analyze the emergent magnetoelectric coupling, demonstrating how it allows us to electrically drive transitions between different topological spin textures employing an interlayer bias.

II. RESULTS AND DISCUSSION
We start our analysis by describing the magnetic order that emerges in twisted CrX 3 bilayers. The magnetic behavior of CrX 3 monolayers can be described by a spin Hamiltonian in a honeycomb lattice (1) where J is the first neighbor intralayer ferromagnetic exchange, taking a value on the order of 2-3 meV [42][43][44] .
A v is the anisotropic magnetic exchange. A u is the singleion anisotropy, which is dominated by A v in CrI 3 , CrBr 3 , and competing with A v in CrCl 3 45 . Our analysis will focus on the anisotropy regime leading to out-of-plane easy axis, that corresponds to CrI 3 , CrBr 3 and slightly strained CrCl 3 4647 . The term V contains additional contributions that do not qualitatively affect our analysis, including Kitaev exchange 48 , biquadratic exchange 49 , direct Dzyaloshinskii-Moriya interaction 50 , and dipolar coupling 51 .
Considering now two layers of CrX 3 stacked together an interlayer magnetic exchange J M will arise. A moiré pattern like the one shown in Fig. 1a emerges for twist angles lower than 3 • . Depending on the stacking between layers two different regions can be distinguished, monoclinic and rhombohedral. Associated with these different regions the sign of J M will change, i.e., J M is ferromagnetic (positive) in the rhombohedric stacking and antiferromagnetic (negative) in the monoclinic one [52][53][54][55] . This leads to a modulation of J M in the moiré supercell like the one shown in Fig. 1b [34][35][36][37]5657 . Therefore, to model the twisted CrX 3 bilayer we add the interlayer interaction to the Hamiltonian of Eq. (1) as where J M (r i , r j ) is the site-dependent interlayer exchange 58 . Since CrX 3 is composed of Cr 3+ with a spin state S=3/2, we can solve in a classical way the spin Hamiltonian for the twisted system 59 . The ground state magnetic order is depicted in Fig. 2a. We can see that a non-collinear magnetic texture emerges between ferromagnetic and antiferromagnetic regions in agreement with previous theoretical 38,39 and experimental results [34][35][36][37]40,41 . Taking this as the starting point, we now show that in the presence of spin-orbit coupling, this topologically-trivial spin texture leads to the emergence of an electric polarization P ij between firstneighbor spins in the same layer separated by a distance r ij due to the inverse Dzyaloshinskii-Moriya (DM) mechanism 60,61 where λ SOC is a coefficient that controls the strength of the spin-orbit coupling and α is a proportionality constant that depends on the electronic structure and crystal environment, which is similar for the three chromium trihalides. The polarization emerges at the middle point between the two neighboring spins S i and S j . From Eq.
(3) we can clearly see the requirement of non-collinearity and the presence of SOC to produce an electric polarization. For the emerging ground state spin texture of the twisted system (Fig. 2a), the associated electric polarization when SOC is introduced is shown in Fig. 2b, with P z the dominant component. From Figs. 2ab, we can observe that an electric polarization emerges in both layers in the areas where the non-collinear spin texture occurs. The polarization emerges locally, and for the ground state spin texture the net electric polarization is zero. Therefore, to analyze the strength of the multiferroic order as a function of the parameters that appear in the spin Hamiltonian, we will consider an average of the electric polarization module P . The module of the electric polarization in the moiré system associated with the ground state spin texture can be seen in Fig. 1c. We can observe there, that the ground state breaks the C 6 symmetry of the original system leading to a C 2 symmetry. The lifted C 6 symmetry stems from a spontaneous symmetry breaking of the ground state due to the ratio between the parameters entering the spin Hamiltonian, analogous to the symmetry breaking associate to stripy magnetic orders. Figure 2c shows the average polarization as a function of the maximum interlayer exchange J max M . At low values, the twisted system behaves like two independent layers, remaining ferromagnetic, no spin texture emerges and consequently, there is no electric polarization. By increasing the interlayer coupling, non-collinearity appears and with it an associated electric polarization in virtue of Eq. (3). Figure 2d shows the average polarization as a function of the anisotropic exchange A v . At high values, the twisted system tends to align the magnetic moments out of the plane. Therefore, no spin texture occurs and thus there is no electric polarization. This situation happens for A v /J > 0.1, so CrI 3 would be on the verge of displaying a multiferroic behavior due to its strong uniaxial anisotropy 45 . We now demonstrate and quantify, using density functional theory calculations 62 , the emergent electric polarization that we have seen that appears in the noncollinear moiré system in virtue of eq. (3). Performing first-principles calculations in a full twisted CrX 3 bilayer with spin-orbit coupling and non-collinear magnetism is well beyond the current computational capabilities. However, since the electric polarization arises locally, the ab initio analysis can be performed in a system like the one shown in Fig. 3a by imposing a magnetic texture in the CrX 3 layer like the one shown in Fig. 3b, that resembles the kind of spin texture found in the ground state between rhombohedric and monoclinic regions. We set the same out-of-plane spin texture to the three compounds to systematically extract the effect of the halide atom and quantification of the inverse DM interaction 63 . The associated polarization of a noncollinear texture given in eq. (3) can be rewritten in terms of a spin spiral propagation vector q and the spin rotation axis e = (0, −1, 0), leading to an electric polarization of the spin texture 60,61 where β is a proportionality constant that depends on the electronic structure and crystal environment and are similar for the three compounds. To demonstrate and quantify the emergence of an electric polarization we performed ab initio calculations in two equivalent magnetic configurations (Fig. 3b) with the same spin rotation vector e, but with opposite spin propagation vector q. Therefore, an opposite electric polarization will emerge in each of the configurations. The emergence of the electric polarization in the spin texture is directly obtained by taking the difference between the two equivalent configurations. This procedure provides a direct methodology to extract electric polarization stemming from non-collinear magnetic order 64 . Modifying the q-vector of the spiral, i.e. the size of the supercell, will change the total value of the ferroelectric polarization obtained since it will modify the non-collinear spin order controlling the inverse DM interaction. In our analysis, we use the same spin texture for every chromium trihalide, and hence we can extract the contribution coming from the spin-orbit coupling to the inverse DM interaction.
The emergence of an electric polarization is accompanied by a reconstruction of the electronic density ρ. Therefore, we can analyze the difference in the electronic density δρ between both configurations. This is shown in Fig. 3c for the three chromium halides. We observe that the electronic reconstruction increases by taking a heavier halide. Thus, for the same spin texture, we can see that CrI 3 will produce the strongest ferroelectric polarization. This result is a consequence of the increase of the spin-orbit coupling when one goes down in the halide group, thus demonstrating the effective equations (3) and (4) governed by the SOC prefactor λ SOC .
The reconstruction of the electronic density will lead to the appearance of a ferroelectric force in the atoms in the spin texture in the direction of the emergent electric polarization 10 . Therefore, we can compute the force difference between both configurations in Fig. 3b to provide a direct quantification of the electric polarization in CrX 3 . Figure 3d shows the force difference between both configurations. We can see that the forces emerge only in the Cr atoms in which the non-collinear magnetism is present. Moreover, we can see that the forces emerge in the z-direction, indicating the direction of the electric dipole, and as expected from the schematic in Fig. 3b. We can quantify the dependence on the halide by taking the average of the force module among Cr-atoms for each of the compounds, as shown in 3d. Taking a heavier halide produces an increase in the ferroelectric force, as expected from the increase of the spin-orbit coupling. As a reference, NiI 2 , a 2D multiferroic governed by this same mechanism of spin-orbit coupling and non-collinear mag-netism, leads to ferroelectric forces of ≈ 20 meV/Å 10 . Therefore, this confirms that spin textures produced in twisted CrX 3 bilayers lead to the emergence of an electric polarization with strong magnetoelectric coupling.
We now elaborate on some conclusions on the strength of multiferroic order in twisted CrX 3 bilayers considering together the results obtained from the low energy model and the ab initio calculations that allowed us to provide a quantification of the inverse DM interaction in these magnetic moiré systems 65 . On the one hand, our ab initio DFT methods show that as the halide becomes heavier, the ferroelectric polarization becomes stronger. On the other hand, the low energy model shows that the anisotropic exchange has a detrimental impact on the formation of the spin texture and hence on the multiferroic order. In particular, in CrI 3 the strong anisotropic exchange might partially quench the formation of a sizable non-collinear magnetic texture. In CrCl 3 , the small SOC would yield a comparably weak ferroelectric polarization despite the formation of the non-collinear texture. Therefore, the optimal twisted bilayer that yields the strongest ferroelectric polarization would be CrBr 3 , or ultimately, a bilayer of an intermediate stoichiometry between CrBr 3 and CrI 3 45 . Finally, we analyze the magnetoelectric coupling present in twisted CrX 3 bilayers 66 . To do so, we now include in the low energy Hamiltonian a coupling to an external electric field E = (0, 0, E z ) perpendicular to the twisted system in the z-direction We now show how this interlayer bias allows controlling the magnetic order due to the emergent multiferroicity 67 . Figures 4a and 4b show the spin texture and the associated electric polarization at E z /J = 3. We can observe that the ground state spin texture shown in Fig.  2a gets significantly modified due to the strong magnetoelectric coupling, leading to the formation of a magnetic skyrmion around the AA rhombohedral stacking in one of the layers 38,39,[68][69][70][71][72][73][74] . Interestingly, such a magnetic state features a non-zero total electric polarization in the z-direction. The evolution of the total electric polarization in the moiré supercell in the z-direction as a function of the external electric field is shown in Fig. 4c. The different bumps in P z as a function of the electric field indicate transitions to different non-trivial spin textures which are shown as insets in Fig. 4c. Starting at E z /J = 0 in the ground state spin texture, at E z /J = ±3 the magnetic skyrmion shown in Fig. 4ab arises. At E z /J = ±5 a second pair of magnetic skyrmions emerge in the other layer around the AB/BA stacking regions 75 . The critical electric field required to produce a transition to the skyrmionic phase is inversely proportional to the strength of the λ SOC . As a reference, J ≈ 3 meV and αλ SOC ≈ 10 −4 e (e electron charge) in twisted CrBr 3 bilayers, implying that such magnetic transitions can be experimentally driven via gating at voltages of 1 − 10 V 67 . This brings twisted CrBr 3 bilayers as a promising platform for the electric control of non-trivial magnetic textures, and ultimately as a platform for magnetoelectric skyrmionics.

III. CONCLUSIONS
To summarize, we have demonstrated that twisted CrX 3 bilayers develop a multiferroic order due to the interplay between the moiré structure, non-collinear magnetic order, and spin-orbit coupling. We have provided a quantification of the strength of the multiferroic order for this family of compounds, finding that, among the stoichiometric chromium trihalides, CrBr 3 is expected to display the strongest multiferroic and magnetoelectric coupling. Furthermore, we have shown that this strong magnetoelectric coupling allows an electrical control of nontrivial magnetic textures using an interlayer bias. Our results put forward a strategy to design a new family of artificial multiferroics with a strong magnetoelectric coupling based on twisted magnetic van der Waals materials.

DATA AVAILABILITY STATEMENT
All data that support the findings of this study are included within the article (and any supplementary files).