Twist engineering of the two-dimensional magnetism in double bilayer chromium triiodide homostructures

Twist engineering, or the alignment of two-dimensional (2D) crystalline layers with desired orientations, has led to tremendous success in modulating the charge degree of freedom in hetero- and homo-structures, in particular, in achieving novel correlated and topological electronic phases in moir\'e electronic crystals. However, although pioneering theoretical efforts have predicted nontrivial magnetism and magnons out of twisting 2D magnets, experimental realization of twist engineering spin degree of freedom remains elusive. Here, we leverage the archetypal 2D Ising magnet chromium triiodide (CrI3) to fabricate twisted double bilayer homostructures with tunable twist angles and demonstrate the successful twist engineering of 2D magnetism in them. Using linear and circular polarization-resolved Raman spectroscopy, we identify magneto-Raman signatures of a new magnetic ground state that is sharply distinct from those in natural bilayer (2L) and four-layer (4L) CrI3. With careful magnetic field and twist angle dependence, we reveal that, for a very small twist angle (~ 0.5 degree), this emergent magnetism can be well-approximated by a weighted linear superposition of those of 2L and 4L CI3 whereas, for a relatively large twist angle (~ 5 degree), it mostly resembles that of isolated 2L CrI3. Remarkably, at an intermediate twist angle (~ 1.1 degree), its magnetism cannot be simply inferred from the 2L and 4L cases, because it lacks sharp spin-flip transitions that are present in 2L and 4L CrI3 and features a dramatic Raman circular dichroism that is absent in natural 2L and 4L ones. Our results demonstrate the possibility of designing and controlling the spin degree of freedom in 2D magnets using twist engineering.

Moiré superlattice forms when two vertically stacked atomic crystals are rotated with respect to each other, enabling a powerful venue to design and tailor the physical properties of 2D materials, including electronic, phononic, and magnetic ones. So far, it has achieved fruitful results in controlling the charge degree of freedom (DoF) and realizing novel quantum phenomena in both single-particle electronic states and two-particle excitonic states. Outstanding examples include the creation of flat electronic bands 1,2 that leads to various strongly correlated 6-10 and topological [11][12][13] phases for the former case, and the introduction of deep confinement potential that localizes excitonic states to realize moiré excitons 14,15 , exciton Mott insulators 16 , and quantum emitters 17 for the latter. Recently, it has also shown significant impacts on modulating the lattice DoF and induced lattice reconstructions 18,19 , renormalized phonons 20 , and moiré phonons 21 . In sharp contrast, the potential and power of moiré superlattices in controlling the spin DoF and engineering the magnetic properties have remained as a pristine area experimentally, despite a few pioneering theoretical predictions including noncollinear magnetism 3 , topological magnetism 4 , moiré magnon bands 5 , and one-dimensional magnons 5 .
Newly discovered layered magnets have greatly expanded the library of two-dimensional (2D) materials and provided exciting possibilities for the exploration and exploitation of the intrinsic spin DoF at the atomically thin limit [22][23][24] . Among the 2D magnetic atomic crystals discovered thus far, atomically thin chromium trihalides, CrX3 (X = Cl, Br, and I) have attracted extensive interest as an archetype 2D magnet platform for hosting a plethora of novel 2D magnetic phenomena [25][26][27][28][29][30][31] and exhibiting versatile tunability with external stimuli [32][33][34][35][36] . In particular, it has been theoretically calculated 37 and experimentally shown 35,36 in bilayer (2L) CrI3 that tuning between the monoclinic (AB') and rhombohedral (AB) stacking geometries can lead to magnetic transitions from the layered antiferromagnetic (AFM) to ferromagnetic (FM) order. Such a close relationship between structural stacking symmetries and interlayer magnetism in CrI3 naturally stimulates the curiosity of investigating the magnetism in CrI3 moiré magnets where the interlayer exchange coupling is periodically modulated in both sign (i.e., AFM vs FM) and strength (i.e., magnitude).
Here, we fabricate twisted double bilayer (tDB) CrI3 homostructures with individual 2L CrI3 to engineer the inter-2L exchange coupling and investigate the resulting moiré magnetism. Our choice of 2L CrI3 as the building block is based on its significantly narrower structural and magnetic phonon linewidths than those of monolayer CrI3, suggesting a much better crystalline and magnetic integrity in 2L CrI3 38,39 . Figure 1a top panel shows a false-color optical image of a typical tDB CrI3 sample that was made by tearing a large piece of 2L CrI3 in the bottom panel into two and then stacking them together at a controlled twist angle of (see Methods). Electron diffraction (Fig. 1b) shows two sets of 1 st and 2 nd order Bragg peaks for the two 2L CrI3 in a tDB CrI3 sample with the targeted of 1.0 o during fabrication, from which the actual was determined to be 0.9 o on average with a standard deviation of 0.1 o through surveying nine different locations on this sample. The match between targeted and measured values confirms our well control of the twist angle in fabrications, and the small standard deviation suggests the homogeneity of our samples. Dark field transmission electron microscopy (TEM) (Fig. 1c) displays the real-space periodic superstructures with noticeable domain formation. The bright triangles are associated with the strongly coupled regions between two 2L CrI3 whereas the dark boundaries represent the decoupled boundaries. Similar reconstructuring has been observed in low-twist angle graphene bilayers 19 and transition metal dichalcogenide bilayers 18 .
Magneto-Raman spectroscopy can capture the interlayer magnetism in few-layer CrI3, by detecting the unique static magnetism-coupled phonons that break the time-reversal symmetry and have antisymmetric Raman tensors [38][39][40][41][42] , in addition to the conventional pure structural phonons. Figure 1d presents representative Raman spectra of tDB CrI3 with the targeted = 1.1 o in both the crossed and parallel linear polarization channels at 10 K, featuring key Raman modes in three spectral ranges, 75 to 85 cm -1 , 95 to 120 cm -1 , and 120 to 133 cm -1 . These modes are coarsely comparable to those of few-layer CrI3 among which the Raman modes in the 75 to 85 cm -1 and 120 to 140 cm -1 ranges particularly highlight the contribution from the magnetism-coupled-phonon scattering 38 . In this work, we focus on the 120 to 133 cm -1 range because the Raman modes here could be related to the moiré magnetism and are of higher intensity than those in 75 to 85 cm -1 . Figure 1e zooms into the Raman spectra in the 120 to 133 cm -1 range and includes data in both linearly crossed and parallel channels with incident polarizations at the horizontal and 45 o rotated directions, i.e., inc = 0 o and 45 o . Clearly, the one primary mode in the parallel channels (U in the triplet in the crossed channel increases in intensity and the two modes on its side (U 2,4 t ) decrease as the twist angle increases, showing the trend that the magneto-Raman spectra of tDB CrI3 evolves from resembling most the natural 4L CrI3 at the lowest to converging towards the 2L CrI3 at the highest .
This trend is expected because the inter-2L coupling strength in tDB CrI3 weakens at larger , leading to the parallel stacked ( = 0 o ) tDB CrI3 relaxing to a 4L CrI3 flake and the large twist angle ones are equivalent to two decoupled 2L CrI3 films. This intensity evolution is further quantitatively summarized in Fig. 2c where the relative intensity ratio U 3 t /( U 2 t + U 4 t ) is plotted against for tDB CrI3 and compared with those of 4L and 2L CrI3, showing a monotonous enhancement with increasing and further confirming the corresponding reduction of inter-2L coupling at larger . Equally informative is the frequency shift of the Raman modes which is shown in Fig. 2b. At low twist angles (e.g., = 0.5 o , 1.1 o ), the frequencies of U 2 t and U 4 t in tDB CrI3 match with those of U 2 4L and U 4 4L in 4L CrI3, respectively, whereas the frequency of U 3 t is close to that of U 2 2L in 2L CrI3 which appears in the crossed channel (the same as U 3 t ) and that of U 3 4L in 4L CrI3 which is however absent in the crossed channel (in contrast to U 3 t ). As the twist angle increases, the frequencies of U 2 t and U 4 t blueshift towards their high-frequency neighbors U 1 t and U 3 t . Eventually, at large twist angles (e.g., = 5 o ), the frequencies of U 2,4 t become nearly indistinguishable from U 1,3 t and approach those of U 1,2 2L in 2L CrI3. The correspondence of mode frequencies between the fabricated tDB CrI3 and the natural 4L/2L CrI3 reveals that magnetism-coupled-phonon contributions for U 2 t and U 4 t arise from regions with strong inter-2L coupling, resembling the 4L-like case, whereas that for U 3 t is dominated by the decoupled regions, mimicking the 2L-like case. The twist angle dependencies of mode frequencies and the relative intensity ratio both confirm the reduction of inter-2L coupling and the suppression of strongly coupled 2L-2L regions in tDB CrI3 at larger .
Out-of-plane magnetic field (B ⊥ ) is known to introduce sharp spin-flip transitions in few-layer CrI3, which is nicely captured in the B ⊥ dependence of Raman spectra of the magnetism-coupled phonon scattering 38 . We now proceed to more in-depth investigations of the engineered magnetism in tDB CrI3 by performing its magnetic field dependencies for selected twist angles and comparing them across one another disappears concurrently at B c = 0.7 T with B c for the spin flip transition from ↑↓ to ↑↑ 25 .
We now turn to the magnetic field dependencies of U of the four Raman modes because U 2 t is overwhelmed by the spectrally closest and much stronger U 1 t in the same co-circularly polarized channels. Clearly, U 1 t and U 3 t show substantial differences between the LL and RR channels with opposite relative intensities, whereas U 4 t is almost helicity independent. Figure 4c presents the B ⊥ dependencies of the fitted intensities of U 1,3,4 t in both co-circularly polarized channels, where the field is swept from + 2 T to -2 T and then back to + 2 T. The two key features highlighted in   46,47 , so that the periodic moiré interlayer exchange interaction is the leading order magnetic energy scale.

Growth of CrI3 single crystals
Single crystals of CrI3 were grown by the chemical vapor transport method. Chromium power (99.99% purity) and iodine flakes (99.999% purity) in a 1:3 molar ratio were put into a silicon tube with a length of 200 mm and an inner diameter of 14 mm. The tube was pumped down to 0.01 Pa and sealed under vacuum, and then placed in a two-zone horizontal tube furnace whose two zones were raised up slowly to 903 K and 823 K for 2 days and then held for another 7 days. Shiny, black, plate-like crystals with lateral dimensions of up to several millimeters were obtained from this growth procedure.

Fabrications of 2D CrI3 and tDB CrI3
Atomically thin 2D CrI3 flakes were exfoliated in a nitrogen-filled glovebox, and their thickness was first determined by the optical color contrast to select natural 4L and 2L CrI3 and further confirmed using Raman spectroscopy at 10 K. Using a polymer-stamping technique inside the glovebox, large-size (lateral dimensions greater than 10 μm) 2L CrI3 flakes were torn into two parts with similar sizes, one of which was rotated by a well-controlled rotation micrometer for targeted twist angles and then brought down to stack with the remaining half. Both 4L/2L and tDB CrI3 samples were sandwiched between two few-layer hBN flakes to avoid surface reactions with oxygen and moisture in the ambient environment after taking out from the glovebox. The samples for magneto-Raman spectroscopy measurements were placed onto the SiO2/Si substrate, and those for TEM measurements were transferred onto to TEM grids.

Supplemental Information for
Twist engineering of the two-dimensional magnetism in double bilayer chromium triiodide homostructures

Background preparations
Because the polar magneto-Raman geometry (i.e., normal incidence and backscattering geometry with the out-of-plane magnetic field) used in the current experiment, one could only access the out-of-plane components of spins due to the optical selection rule. Therefore, in our model analysis below, we will only focus on the out-of-plane components of spins. In principle, the frustration among spins in moiré magnets could result in spin canting off the out-of-plane direction even in Ising-type magnets 1 . While we do not explicitly include possible in-plane components in our model, it is worthwhile pointing out that the main conclusions from the model calculations below are not sensitive to the presence or absence of spin canting off the out-of-plane direction and are compatible with the noncollinear magnetism.
In a moiré superlattice, the periodic modulations of the interlayer stacking geometry lead to the magnetic domain structures in twisted CrI3 2 . While the structural stacking pattern and the resulting magnetic domain structure are sophisticated, it has been shown theoretically that two types of stacking geometries, AB (i.e., rhombohedral) and AB' (monoclinic) stacking, have the lowest elastic energy 3 , and therefore, shall be the dominant ones. Instead of trying to capture all microscopic details about the stacking pattern and the lattice reconstruction, in our theory model, we will take a simplified setup and only focus on the most relevant and dominant parts: the monoclinic AB' and the rhombohedral AB stacking, while omitting the other highenergy stacking configurations. In reality, these high-energy stacking area may arise and fill up the regions other than the monoclinic and rhombohedral areas, especially at twist angles greater than ~ 1.0 o 2 . But, the omission of them is sufficient to the leading order approximation.
In the twisted double bilayer (tDB) CrI3 samples studied in our current work, the interlayer exchange coupling within individual bilayer (2L) CrI3 remains antiferromagnetic (AFM), just the same as that in regular isolated 2L CrI3, for which we take the value from the literature, 1 = 0.04 meV/μ B 2 2,3 . The coupling between the two 2L CrI3 (i.e., at the interface between the two 2L CrI3) requires the consideration of two types of structural stacking: monoclinic AB' stacking that corresponds to interlayer AFM exchange coupling and rhombohedral AB stacking that has interlayer ferromagnetic (FM) exchange coupling. We again adopt the values from literature 2,3 , 1 = 0.04 meV/μ B 2 for AFM and 2 = −0.6 meV/μ B 2 for FM. Please note that, according to the first-principles calculations 2,3 , | 2 | is significantly greater than | 1 |. As it will be shown later, this fact is crucial for the development of net magnetization in the tDB CrI3 with intermediate twist angles. Finally, the intralayer exchange coupling is known to be FM ( < 0) and the firstprinciples calculations give = −2.2 meV/μ B 2 .

Theoretical modeling in tDB CrI3
For the monoclinic inter-bilayer stacked regions, the interlayer coupling between any two neighboring layers are all AFM with 1 = 0.04 meV/μ B 2 , and therefore its magnetism favors the interlayer AFM order: ↑↓↑↓ or ↓↑↓↑ from the top to the bottom layer. In contrast, for the rhombohedral inter-bilayer stacked regions, the preferred spin alignment is ↑↓↓↑ or ↓↑↑↓ because of the FM exchange coupling at the bilayerbilayer interface. It is important to note that for both stacking geometries, all these energetically favored spin configurations have zero net magnetization. In the other word, neither of the stacking favors the development of a net magnetization that is observed in our experiment: the large Raman circular dichroism for the 1.1 o twisted tDB CrI3 at 0 T. As we shall see below, to achieve the experimentally observed net magnetization in 1.1 o tDB CrI3 at 0T, it requires the frustrated spin interactions due to the competition of spins between the monoclinic AFM and the rhombohedral FM regions.
 the magnetic state with a zero total magnetization Case I: One very obvious possible configuration is that the monoclinic region has ↑↓↑↓ (or the opposite one), whereas the rhombohedral domain has ↑↓↓↑ (or the opposite one). In this configuration, each stacking region has its lowest-energy configuration, but the boundaries between them have opposite spin alignments within each of the bottom two layers. The energy penalty associated with such (intralayer) magnetic domain walls is proportional to the length of the domain boundary L multiplied by a factor of 2 because of their presence in two layers and the intralayer FM exchange coupling , i.e., 1 ∝ 2 / where is the lattice constant of CrI3. Because the intralayer coupling strength is much greater than the interlayer ones, this energy penalty will become large and dominant, when the domain boundary becomes long.
Case II: If the energy cost for such intralayer domain boundaries becomes too high, the system will favor a different configuration, i.e., avoiding any intralayer magnetic domain boundaries. Take an extreme case example, if the area of the monoclinic domain is much larger than the rhombohedral one, we would expect ↑↓↑↓ for both monoclinic and rhombohedral stacking regions. This spin configuration does not have any intralayer domain wall penalty, although the rhombohedral one is not in its ground state and is subject to the energy penalty of 2 ∝ | 2 | ( / ) 2 , where is the length of the domain boundary and the area shall scale with 2 .
Depending on the size of the domain , which is controlled by the twisting angle of tDB CrI3, the relative ratio between 1 and 2 shall change correspondingly. Most importantly, when the domain size is large (small), we shall expect 1 ≪ 2 ( 1 ≫ 2 ), and thus the spin configuration shall follow the first (second) possibilities described above. Here, we emphasize again that all these spin configurations have zero total magnetization.

 the magnetic state with a non-zero net magnetization
Case III: For the intermediate domain size, where 1~2 , new possible spin configuration can arise, for example, ↑↓↑↓ for the monoclinic region and ↑↑↑↓ for the rhombohedral region. In this case, the rhombohedral region is not in its ground state, because the top two layers (i.e., the top 2L CrI3) both have spin up that is in contrast with the interlayer AFM coupling for a 2L CrI3. This results in an energy penalty of ~1 ( / ) 2 . At the same time, the intralayer magnetic domain boundary only exists within one layer (i.e., the 2 nd layer), which costs an energy of ~/ . Therefore, in total, the energy cost for such a spin configuration is 3 ∝ 1 ( / ) 2 + / . Please note that this configuration has only one layer with intralayer magnetic domain wall (in the 2 nd layer), so its intralayer energy cost is less than the Case I spin configuration where two layers are subject to intralayer magnetic domain wall (in the 3 rd and 4 th layers). This is the reason why this Case III configuration is possible to be more energetically favored than the Case I spin arrangement. In comparison to the Case II spin configuration, this Case III configuration has smaller interlayer energy cost, because 1 < | 2 |.
This Case III spin configuration contains both intralayer and interlayer energy penalty ( 3 ), but the intralayer part is smaller than that of Case I ( 1 ) and the interlayer one is smaller than that of Case II ( 2 ). Therefore, when 1~2 , this configuration could be most energetically favored.

 the comparison of energies for the three spin configurations
Below we compare the energies of these three configurations (Case I, II, and III) as a function of the domain size . Here, we assume that the monoclinic region is greater, because it is the natural stacking for fewlayer CrI3. Here, we use to represent the linear size of the rhombohedral stacking region (normalized by the lattice constant ), and 2 to represent its corresponding area. In reality, it should be multiplied by a geometric factor of order 1, which depends on the shape and other geometrical details of the moiré superlattices, and this factor is set to 1 here for simplicity. At a qualitative level, our conclusion is robust to the exact value of this factor. For a similar argument, we set the length of the magnetic domain boundary to be 4 .
From Fig. S1, we can see that there are two critical lengths, L1 (~ 4 / 1 ) and L2 as marked, at which the most energetically favored spin configuration alters. Let's discussion the three regions below, above L1, between L1 and L2, and below L2.
(a) Above L1 when the rhombohedral region dominates over the monoclinic region (i.e., not quite practical -the parameter range where this simple model fails when the twist angle is so big that the inter-bilayer coupling becomes so weak), the most energetically favored spin configuration is Case I: the monoclinic region has ↑↓↑↓ (or the opposite one) and the rhombohedral domain has ↑↓↓↑ (or the opposite one). This gives zero total net magnetization. (b) Between L1 and L2 when the rhombohedral and monoclinic regions are comparable (i.e., intermediate twist angle), the most energy-saving spin configuration is Case III: the monoclinic region has ↑↓↑↓ and the rhombohedral domain has ↑↑↑↓. This configuration gives a finite total net magnetization, which is consistent with our result in the 1.1 o tDB CrI3. (c) Below L2 when the monoclinic region dominates over the rhombohedral region (i.e., very small twist angle), the most energetically favored spin configuration is Case II: homogeneous across the monoclinic and rhombohedral regions ↑↓↑↓. This again gives zero total net magnetization.

Figure S1
Plots of calculated energy penalty ( 1 , 2 , and 3 ) for the three spin configurations as a function of the rhombohedral stacking region size, normalized to the lattice constant a, i.e., L/a.

 the consistency between calculations and experiment for the 1.1 o tDB CrI3
We can do a quantitative comparison between our theoretical calculations and experimental observation. It should be emphasized due to the complicated nature of this material, we made approximations and simplifications in the theory treatment above. In addition, our theory has no fitting parameter at all, and the values of all control parameters are obtained from first-principles calculations. On the other hand, for small energy scales, like interlayer spin exchange, first-principles techniques are limited by resolution and accuracy, which is another source of error for our predicted values. Hence, the theory predicted values here should be treated as an order of magnitude estimation, instead of an exact solution.
With these limitations in mind, here we give key values from the theory predictions. The Case III spin configuration with a net magnetization shall arise at the critical domain size of 1 ≈ 4 / 1 . Using values from first-principles calculations 2,3 , ≈ 220. This value is expected to scale with the linear size of a moiré unit cell. For our intermediate angle at which we observed Case III in the experiment, = 1.1°, we expect ≈ 1/ ≈ 52, which satisfies the requirement of < 1 and on the right order of magnitude with 1 .
Moreover, because the rhombohedral region has net magnetization, under an external magnetic field along (against) the direction of this net magnetization, the magnetic domain for the rhombohedral stacking shall increase (shrink) in size to reduce the energy. This magnetic domain size change is a smooth function of the magnetic field strength, and this is the reason why the signal intensity observed in this phase can vary continuously, gradually as the magnetic field increases.

S2. Calculations of the magnetism-coupled phonon contributions in tDB CrI3
In few-layer CrI3, it has been established that the interlayer structural coupling can lead to Davydov splitting of the Ag phonon mode of CrI3, and these split phonon multiplets have two contributions in the Raman scattering: the pure structural and layered magnetism-coupled phonon contributions that correspond to the fully symmetric and antisymmetric Raman tensors 4 . Here, the layered magnetism-coupled phonon contribution encodes the static magnetic order information and tracks its evolution as a function of external magnetic field.
Here, we exploit the same idea for the strongly coupled 2L-2L regions, including both the monoclinic and rhombohedral inter-2L stacking regions, in tDB CrI3. We perform the same analysis as in Ref. 4 to compute the magnetism-coupled phonon contributions for the various magnetic orders in the two stacking regions in different magnetic field ranges. The results are summarized in Table. S1 below.
Please note that the atomic displacement eigenvectors are nearly identical for the monoclinic and rhombohedral stacking regions, sharing the same symmetry properties although having slightly different atomic displacement magnitudes across the layers. Therefore, they are listed with the same sketch in Table. S1. Table S1. Magnetism-phonon coupling strength ( • / ) between the phonon mode ( ) and the magnetic order ( / ) in either monoclinic (M in the superscript, red) or rhombohedral (R in the superscript, blue) stacking regions in tDB CrI3, across all four magnetic field regions with different out-ofplane spin configurations. In each entry,√ (×) stands for the presence (absence) of the magnetism-phonon coupling.

Modes
Magnetic States