Magnetization textures in twisted bilayer 2D CrX$_3$ (X=Br, I)

Motivated by the recent experiment demonstration of stacking dependent interlayer magnetic interaction [T. Song et al., Nat. Mater. 18, 1298 (2019); T. Li et al., Nat. Mater. 18, 1303 (2019); W. Chen et al., Science 366, 983 (2019)], we investigate the magnetization textures and the control possibilities in the moir\'{e} pattern formed of twisted bilayer two-dimensional (2D) magnets CrX$_3$ (X=Br, I). We find that the stacking dependent interlayer magnetic interaction results in the formation of periodic magnetization domains in a long-period moir\'{e} pattern. Magnetization textures with various topological numbers can be constructed, depending on the winding of the textures around the domain walls. A uniform external magnetic field competes with the lateral modulated interlayer magnetic interaction and can be utilized to tune the magnetization textures.

The physical properties of vdW layer materials depend sensitively on the interlayer interaction, which is further connected to the interlayer stacking order of the atomic structure. For example, in bilayer graphene, the lowenergy state in AA configuration features linear dispersion, while the one in AB configuration is parabolic [36]. In bilayer transition mental dichalcogenides, both resonance Raman and photoluminescence spectra show distinct features in 3R and 2H stacking configurations [37]. The stacking dependence is more prominent in vdW 2D magnets, in which the stacking order is directly related to the arrangement of magnetic momentums of the two layers. Indeed, in bilayer CrI 3 , first-principles studies have revealed a rich phase diagram between interlayer * Electronic address: tongqj@hnu.edu.cn ferromagnetic (FM) and antiferromagnetic (AFM) orders as the stacking configuration is continuously changed [12,[38][39][40][41]. This interesting stacking dependent interlayer magnetism has been confirmed experimentally in CrI 3 [42,43], and also demonstrated in CrBr 3 [44], suggesting a way for designing magnetism via controlling stacking order. The stacking dependent interlayer magnetic coupling is expected to be general in magnetic vdW materials, as the magnetic interaction depends intimately on both relative distance and orientation between the magnetic moments.
Moiré pattern ubiquitously forms in vdW layered materials due to the misorientation or lattice mismatch between the constituent layers [45][46][47][48][49][50][51][52][53][54][55]. In a long-period moiré pattern, the interlayer atomic configuration resembles lattice commensurate structure, while the stacking order changes smoothly over long range. Because the interlayer interaction generally depends on the stacking order, such a spatial change of the stacking order then introduces a lateral modulation of interlayer interaction in the moiré pattern. It has been suggested that such a spatially changing interlayer magnetic interaction in 2D magnets can be used to generate various magnetization textures [56][57][58][59]. Extending these ideas to realistic material system would facilitate their experimental observation and potential application in device fabrication.
In this work, we show that, nonuniform magnetization textures can be formed in the moiré pattern of twisted bilayer 2D ferromagnet CrI 3 and CrBr 3 arising from the stacking dependent interlayer magnetic interaction. In a moiré unit cell, there are three magnetization domains surrounded by the magnetization textures of opposite direction. These magnetization textures can be engineered to have topological number ranging from 0 to 3. An external magnetic field competes with the stacking dependent interlayer magnetic interaction and can switch the three magnetization domains to the ones of opposite direction, creating a symmetric hysteresis with double loops in the magnetic dynamics. Arising from the topological protection, the topological number is conserved during the whole process. These interesting magnetization textures can be generated via annealing from a paramagnetic state. Our results suggest new possibili- ties for the investigation of 2D magnetism and suggest a route towards nanoscale magnetization textures by moiré engineering.
The rest of the paper is organized as follows. In Sec. II, we give a brief account of atomic structure of the moiré pattern in bilayer ferromagnets. The stacking dependent interlayer magnetic interaction is also given by first-principles calculations. Different magnetization textures and their topological properties are present in Sec. III. The magnetic field dependence of the magnetization textures and the dynamics are given in Sec. IV. Finally, a discussion on the preparation of the magnetization textures and a summary are given in Sec. V.

II. ATOMIC STRUCTURE OF THE MOIRÉ PATTERN AND STACKING DEPENDENT INTERLAYER MAGNETIC INTERACTION
Monolayer CrX 3 (X=Br, I) has a honeycomb lattice structure possessing a magnetic moment of 1.5µ B per Cr atom [60]. For bilayer CrX 3 , a small twisting or strain between the layers creates a long-period moiré pattern, as schematically shown in the Fig. 1. The moiré periodicity is approximately A ≈ a/ √ δ 2 + θ 2 for small lattice mismatch δ and/or twisting angle θ, where a is the lattice constant of the monolayer. Because the periodicity A can be further tuned by a relative twisting or strain between the layers, it is taken as a variable in the following. In a long-period moiré pattern, the stacking order in each local region R is similar to the lattice-matched stacking configuration but changes smoothly over long range. Typical bilayer configurations are shown in Fig. 1, which are named as AB 1 , AC , AB , AB, AA, and BA stacking in the following. The AB and BA stackings are the same as rhombohedral stacking, and the AB -stacking is the same as monoclinic stacking [61]. In Ref. [44], a point close to AB 1 stacking in CrBr 3 is named as R-type stacking. The interlayer stacking order at any local region R in the moiré pattern can be characterized by an interlayer translation vector r, which is defined in a monolayer unit cell (c.f. Fig. 1(b)). The stacking configurations of AB 1 , AC , AB , AB and AA correspond to the top layer laterally shifted by r = a1+a2 with respect to bottom layer originated from BA configuration, where a 1 and a 2 are the two unit vectors of the monolayer. In a moiré pattern, the number of atoms involved is ∝ (A/a) 2 , therefore a direct calculation of its magnetic property is challenging. However, because the atomic registry in the moiré pattern changes smoothly, each local region resembles lattice matched commensurate structure. Based on this observation, in the following we first study the interlayer magnetic interaction at each local region and then extend to the whole moiré pattern by introducing an effective magnetic field.
The interlayer stacking-order dependent magnetic interaction is related to energy difference between interlayer FM and AFM spin configurations, which can be obtained from first-principles calculations and are shown in Figs. 2(a) and 2(b) for CrI 3 and CrBr 3 respectively. The details of the first-principles calculations and structure parameters are given in Appendix A. One can see that the ground state modulates between FM and AFM as the interlayer stacking order changes. For both CrI 3 and CrBr 3 , the three C 3 -rotation symmetric stackings (AA, AB and BA) favor FM state. However, the magnetic interlayer interaction for C 3 -rotation symmetry broken configurations are different for these two materials. For CrI 3 , the AB 1 , AB and AC regions favor AFM ground state, while for CrBr 3 , the AB 1 region favor AFM ground state, however AB and AC regions favor FM ground state. These results are consistent well with the recent experiment observation of stacking-dependent interlayer magnetism, where the monoclinic (AB ) stacking in CrI 3 [42,43] and R-type (AB 1 ) stacking in CrBr 3 [44] has been confirmed to possess AFM ground state. Figs. 2(c) and 2(d) give a detailed variation of the energy difference between interlayer FM and AFM states along the [100] (blue-dashed line with dots) and [110] (red-dashed line with squares) direction of a unit cell, which includes all the above six stacking configurations. One can see that the magnitude of interlayer magnetic interaction of CrI 3 is stronger than that of CrBr 3 . Such a stacking order dependent interlayer magnetic interaction, combined with the dramatic tunability of the former by a mechanical means, point to exciting possibilities to tailor the magnetic states in the atomically thin vdW 2D magnets [42][43][44].

III. MAGNETIZATION TEXTURES IN THE MOIRÉ PATTERN
With the knowledge of stacking dependence of interlayer magnetic interaction in the commensurate bilayers, we now turn to study the possible magnetization textures in the moiré pattern of 2D magnets. We assume that the magnetic order in one of the two layers is fixed, for example by a substrate. We name the two layers as fixed and free layers in the following. Then at different locals the magnetic order in the free layer tends to point at different directions depending on the interlayer magnetic interaction. Explicitly, if the interlayer interaction is FM, the magnetic order in the free layer tends to align parallel with the one in the fixed layer. On the other hand, if the interlayer interaction is AFM, it tends to align antiparallel with the one in the fixed layer. The effect of a fixed layer on the free layer then can be regarded as an effective magnetic field B(r), which is stacking order r dependent. This effective field can be calculated from where E FM and E AFM are the energies of the interlayer FM and AFM states in a unit cell and m = 1.5µ B z is magnetic moment of Cr atom. Note that the effective field is stacking dependent, which aligns with the magnetization of the fixed layer when the ground state favors FM order, while anti-aligns when the ground state favors AFM order.
In a moiré pattern, the interlayer stacking order changes smoothly, and one expects that the spatially modulated interlayer magnetic interaction (c.f. Fig. 2) then defines a spatially modulated effective magnetic field for the free layer. We assume that the moiré pattern is formed by two rigid lattices (i.e. no lattice reconstruction). In this case, the mapping between the local registry (characterized by the interlayer translation vector r defined in Fig. 1b) and the location R in the moiré pattern is a linear transformation [48,62]. For small twisting angle θ and biaxial strain δ, the mapping functions take We have checked that these two types of mapping functions give almost the same results. In the following, we use the mapping form of biaxial strain to define the moiré magnetic field B(R(r)). Under this local approximation, the magnetization of free layer can be simulated with the effective Hamiltonian including the spatially changing effective magnetic field B(R i ) and intralayer magnetic interactions, where i, j covers all nearest neighboring sites of the hexagonal lattice and m i is magnetic moment at i -th Cr site. J is the intralayer exchange coupling and K is the magnetic anisotropic energy, which are given in the Table II in Appendix A. We have also added a uniform external magnetic field B ext to tune the magnetic configuration. The dipolar interaction is orders weaker than the exchange interaction and neglected here. Via introducing this moiré magnetic field B(R), the bilayer moiré system is simplified to a monolayer honeycomb lattice experiencing a spatially changing external magnetic field of moiré periodicity A, which is discretized on the monolayer honeycomb lattice in moiré scale. The spatially modulated effective magnetic field tends to create nonuniform magnetization domains. While the intralayer exchange and anisotropic interaction favor uniform magnetization texture. Their competition determines the final magnetization configuration. The steadystate and dynamics of the magnetization textures can be solved from the Landau-Lifshitz-Gilbert equation, where H ef f i = − ∂H ∂mi , and γ and α are the gyromagnetic ratio and Gilbert damping coefficient respectively. We obtain various magnetization textures by relaxing from different initial magnetization configurations. The initial magnetization configurations with different topological numbers can be designed via shaping the in-plane magnetic moments around the domain walls to align parallel or form a vortex. α = 1 is used for faster convergence When atomic lattice relaxation is considered, the mapping between the local atomic registry and moiré superlattice would become nonlinear. Furthermore, when the in-plane vector in local registry space becomes divergent, the atomic configuration in moiré space would have rotation [62,63]. We also note that the moiré pattern is quasi-periodic, in which the atomic structures are slightly different in different moiré unit cells. These effects would change the distribution of effective magnetic field B(r) in the moiré pattern. The lattice reconstruction would also change the relative positions of neighboring magnetic moments used for the numerical calculations, which are fixed for the moiré pattern formed by two rigid lattices. Including these lattice reconstruction effects should take into account the competition between the adhesion and the elastic properties of the 2D magnets, which would be an interesting topic for future studies. When the twisting angle is large, the atomic reconstruction is weak and the linear mapping is valid [62,64].
In the moiré of bilayer magnets, we indeed find the formation of nonuniform magnetization textures, including topologically trivial magnetic bubble (MB) and magnetic skyrmion (SK) with different topological numbers, as shown in Figs. 3(c-f). We name the magnetization aligned or anti-aligned with the one in fixed layer as aligned or anti-aligned magnetization in the following. One can see that there are three anti-aligned (red) domains, centered at the AB 1 stacking regions. For CrBr 3 , the uniform FM state is unstable and energetically higher than the MB state for a moiré periodicity larger than 42a, corresponding to a twisting angle of ∼ 1.4 • . For CrI 3 , MB texture forms for a smaller moiré with a periodicity of 18a, arising from the much larger interlayer magnetic interaction (c.f. Fig. 2). Figs. 3(a) and 3(b) show that, as the periodicity increases, the energies of all of the four nonuniform textures decrease. This is because in a long-period moiré, the magnetization orders smoothly that reducing the intralayer exchange energy and is hence dominant by the interlayer moiré magnetic interaction, which has a sign change in the moiré supercell. Therefore, the periodic magnetization domains are more stable in a long-period moiré.
Figs. 3(c-f) present four typical low-energy magnetization textures, which possess three anti-aligned magnetization domains surrounded by the aligned magnetization texture. Around the domain walls, the magnetization textures have different winding manner, i.e. either aligning parallel or forming a vortex. We note that the formation of nonuniform in-plane magnetization in SK state costs more intralayer exchange energy, leading to a higher energy than MB state and the energy increases linearly with the increase of SK number, as shown in Figs. 3(a) and 3(b). Therefore the topologically trivial MB state is always the lowest energy magnetization texture. These nontrivial magnetization textures can generate an emergent electromagnetic field defined as [65] Ω (r) = 1 2 m · ∂m ∂x × ∂m ∂y (4) The lower panels in Figs. 3(c-f) show the emergent electromagnetic field Ω (r), which is mainly distributed around the domain walls manifesting the changing of magnetization textures. One can see that Ω (r) in a domain wall with vortex structure is much larger than the one when parallelly aligned, demonstrating that the former possesses a heavier deformation and costing more intralayer exchange energy. When transporting in such a magnetization texture, carriers would experience a topological Hall effect [65]. The integration over the moiré unit cell gives a quantized topological charge C = 1 2π Ω (r) d 2 r. Our calculation shows that the corresponding topological charges are 0, 1, 2 and 3 for configurations (c-f) respectively. We name these configurations as 3MB, 1SK+2MB, 2SK+1MB and 3SK accordingly. Note that there are other configurations of different topological charge, which are energetically higher and can be found via designing the initial configurations with higher winding number around the domain walls. These higher topological number skyrmions can also be created via annealing from a paramagnetic state with randomly aligned magnetic moments, as discussed in Sec. V.

IV. MAGNETIC CONTROL ON THE MAGNETIZATION TEXTURE
An external magnetic field B ext along z -direction prefers uniform spin polarized texture, which would compete with the spatially changing effective field and provide an efficient means to tune the nonuniform magnetization textures. Fig. 4 shows the steady-state magnetization structure when an external magnetic field is applied for 3SK magnetization texture (see Appendix B for magnetic control of other magnetization textures). With the increase of the magnetic field, the three anti-aligned domains first become larger and then merge into each other, with three aligned domains appearing centered at AA, AB and BA stacking regions for B ext = 0.3J. When further increasing the magnetic field, the aligned domain around AA stacking disappears first, and then the domains at AB and BA regions disappear. This is because the interlayer FM magnetic interactions in AB and BA regions are larger than the one in AA region. For B ext > 1.02J, the magnetization texture becomes a spin polarized state. With the increase of the magnetic field, the energy of magnetization texture (black line with diamonds in Fig. 4) first increases. This is because the Zeeman energy −B ext M z is positive for a negative net magnetization M z (purple line with squares). When M z become positive for B ext > 0.24J, the Zeeman energy become negative and hence the energy deceases with B ext now.
To study in detail the transition from the three antialigned domains to the three aligned domains in the presence of an external magnetic field, we plot the magnetic dynamics when slowly sweeping up (purple-solid line) and then down (orange-dashed line) the magnetic field over a small range (from 0.18J to 0.27J ) in Fig. 5. With an applied uniform external magnetic field B ext , the total effective field felt by the monolayer magnet changes to B ext + B(R), which drives slowly the magnetization dynamics when B ext changes adiabatically. The intralayer exchange and anisotropic energy tend to keep the initial shape of the originally aligned magnetizations. On the other hand, the slowly changing effective magnetic field tends to drive the magnetization to follow its dynamics, which would flip the aligned magnetizations and merge the three disconnected anti-aligned domains. For small external magnetic field B ext , the former effect dominates and the three anti-aligned domains are still disconnected. As the magnetic field increases beyond some critical value, the aligned magnetizations at the merging areas flip and the three anti-aligned domains merge quickly with the formation of disconnected aligned domains. In the reversal process, the anti-aligned magnetizations also tend to keep their original shapes. When the magnetic field decreases beyond some critical value, the anti-aligned magnetizations at the merging areas flip and the aligned domains merge quickly with the formation of disconnected anti-aligned domains. The averaged magnetization M z in the whole process forms a magnetic hysteresis loop. Interestingly, the averaged magnetization M z shown in Fig. 5 forms a double hysteresis loop, with the lower loop ranged at [0.218J, 0.232J ] and the upper loop ranged at [0.232J, 0.245J ]. From the evolution of magnetization texture, one can see that the lower loop forms when parts of the three anti-aligned domains become merged. As B ext increases to the value of 0.232J, the rest of the three anti-aligned domains start to merge, giving rising to the upper loop. When all of the anti-aligned domains become connected, the three aligned domains appear. Different from the vortex domain wall structure of the initial magnetization texture, the ones in the aligned domains are surrounded by one vortex domain wall with winding number of 1 centered at AA stacking and two anti-vortex domain walls with winding number of −2 centered at AB and BA stackings. The main reason of this arrangement of winding number is that, in order to minimize the intralayer exchange energy, the in-plane magnetizations tend to align parallel. With the increase of external magnetic field, the three anti-aligned domains come close and finally merge into each other. In order to keep parallelly aligned at the merging areas, the in-plane magnetizations form a vortex structure with winding number of 1 around the AA stacking domain. While, around the AB and BA domains, the winding numbers become -2. The in-plane magnetizations around the domain wall of AA stacking domain would also tend to align parallel with the nearby in-plane magnetizations around the domain walls of AB and BA stacking domains. This leads to a π/2 helicity of the in-plane magnetizations around the domain wall of AA stacking domain, which is different from the zero helicity of the initial configuration. For initial states with 1SK and 2SK shown in Figs. 3(d) and (e), we have checked that the winding numbers around the AA AB, and BA domains are {0,-1,0} and {0,-1,-1} respectively at the end of the magnetic dynamics. The change of anti-aligned to aligned domains, together with the change of domain wall structures, ensures that the topological charge is preserved during the evolution. When decreasing the magnetic field, the magnetization texture returns to its initial configuration. The whole process then realizes a topologically protected magnetic switching. We note that if further increasing the magnetic field to some critical value, where a topological phase transition happens, the topological charge in the whole process would not be conserved anymore.
For CrBr 3 , which has much weaker anisotropic energy (see Table II) and smaller magnitude of interlayer coupling (see Fig.2 (b)) than CrI 3 , the magnetization domains are easier to break their original shapes, which results in a much smooth hysteresis loop, as shown in Fig.  10 in Appendix C. Because the aligned domains centered at the AB and BA stackings are formed almost simultaneously, there is only one hysteresis loop in this case. Furthermore, the hysteresis loop is easier to observe in a long-period moiré pattern with small twisting angle. For a short-period moiré pattern with large twisting angle, the domain size is small and it costs less intralayer exchange energy for it to merge into each other. Hence, the magnetic dynamics in a short-period moiré pattern is almost reversible, as shown in Fig. 11 in Appendix C.

V. DISCUSSIONS AND SUMMARY
The ground state magnetization texture in a longperiod moiré magnet is MB state, which can be spontaneously prepared by first polarizing the magnetization in a spin polarized state with a strong external magnetic field and then ramping down the field for the magnetization to relax to the ground state. To prepare the SK texture, which is excited state, we first prepare the magnetic state to the paramagnetic state, in which the magnetization aligns randomly and then reduce the temperature to a ferromagnetic state. From this annealing process, one can prepare various magnetization configurations. Fig.  6 shows the simulation results evolving from a random magnetization configuration. Besides the SK with vortex texture of C = 1, there also exists SK with antivortex texture of C = −1, and magnetization textures of high topological number of C = ±2 in each local domain. All of these textures are metastable configurations with energies higher than the MB state. Because the initial paramagnetic state, in which all of the magnetic moments align randomly, has no periodicity, the local topological charges after annealing in the steady state have no local periodicity as well. In a long-period moiré pattern, the local magnetic domains are separate in moiré space and almost independent. Therefore, it is possible to have neighboring magnetization domains with aperiodic topological index. They would animate into topological trivial MB states when domains with aperiodic topological index come close to each other, for example by applying an external magnetic field.
In summary, we have studied the magnetization textures in the moiré pattern of twisted bilayer CrX 3 (X=Br, I). We show that periodic magnetization domains can be formed arising from the lateral modulation of interlayer magnetic interaction. Magnetization textures with various topological numbers can be generated via annealing from a paramagnetic state. An external magnetic field can be used to tune these magnetic configurations. These intriguing magnetization textures with high controllability not only point to new possibilities to study magnetism in 2D limit, but also provide a potential platform to construct scalable spintronic devices.  [66,67]. In order to account for strong electronic correlations for the Cr atoms, a Hubbard onsite Coulomb parameter of 3eV was selected in the calculations [68]. A vacuum layer with thickness of 15Å is used to eliminate the interaction between the layers. The convergence criteria for energy and force are set as 10 −5 eV/Å and 10 −3 eV/Å, the Monkhorst-Pack mesh is set to 12 × 12 × 1, and a plane-wave cutoff energy of 450eV was used in the calculations. The calculated structure parameters of AB-stacking bilayer CrX 3 (X=Br, I) are showed in the Table I. For simplicity, in studying the translation dependence of interlayer magnetic interaction in Fig. 2, the structure parameters are fixed as the one in AB-stacking. The intralayer exchange coupling J and magnetic anisotropic energy K of the monolayer CrX 3 (X=Br, I) are given in Table II. To see how the interlayer magnetic interaction depend on the vdW corrections, we perform calculations using optB88-vdW and optB86-vdW functionals [70]. From     in a long-period moiré pattern.

APPENDIX B. MAGNETIC FIELD CONTROL OF 3MB MAGNETIZATION TEXTURES
In the main text, we have studied the magnetic control of 3SK magnetization texture. We find that other textures show similar magnetic field dependent behavior. In Figs. 8 and 9, we present results for 3MB magnetization texture. Although the topological number is differ-ent, the energy, averaged magnetization and its dynamics show similar behavior as the 3SK magnetization texture. In Figs. 10 and 11, we study the magnetic field control of twisted bilayer CrBr 3 and short-period moiré pattern. The much weaker anisotropic energy and smaller magni-tude of interlayer coupling in CrBr 3 make the magnetization domain easier to break its original shape, resulting in a much smooth hysteresis loop (Fig. 10). For a shortperiod moiré pattern, the domain size is small and it costs less intralayer exchange energy to merge into each other. Therefore the magnetic dynamics in a short-period moiré pattern is almost reversible (Fig. 11).