Twist-assisted all-antiferromagnetic tunnel junction in the atomic limit

Device fabrication

We prepared electrical contacts on Si/SiO₂ (285 nm) substrates by standard electron-beam lithography (Raith PIONEER Two) followed by Ti (2 nm)/Au (10 nm) or Ti (2 nm)/Au (30 nm) deposition (scia Coat 200) to form the electrodes. These electrodes are distributed in a ‘T’ shape. All flakes (HQ Graphene) were exfoliated onto Si/SiO₂ (285 nm) substrates using Scotch Magic Tape. The exfoliation for thin graphite and hBN was carried out under ambient conditions and for CrSBr was carried out under inert conditions in a N₂-filled glovebox with <3 ppm O₂ and <1 ppm H₂O content. The thickness of the CrSBr flakes was initially identified by optical contrast, further examined by Raman spectroscopy (Horiba LabRAM HR Evolution) and atomic force microscopy (Bruker Dimension Icon); see Supplementary Fig. 1. The layer-dependent magnetism deriving from magneto-transport measurements confirmed the thickness of the flakes. The vdW assembly of the hBN/thin graphite/CrSBr/thin graphite/hBN heterostructure was carried out in the same glovebox by sequentially picking up each flake using a poly(bisphenol A carbonate)/polydimethylsiloxane stamp⁴⁴. To insert a hBN monolayer into the twisted interface, an extra pick-up process is required. We used two strategies to align the twist angles. For large exfoliated CrSBr flakes, the tear-and-stack method was used^45,46 (Supplementary Fig. 2a–c). On the other hand, we noticed that the exfoliated CrSBr flakes are always in a ribbon shape with a long edge along the a axis and a short edge along the b axis^33,37, so we aligned small flakes with respect to their long edges (Supplementary Fig. 2d–f). The twist angles based on both of these strategies match well with the experimental results, confirming their reliability. The finished assembly on the stamp was released onto the SiO₂/Si substrate with pre-prepared electrodes and letting the thin graphite electrodes make contact with two of them. The residual poly(bisphenol A carbonate) was dissolved in chloroform before the electrical-transport measurements were carried out.

Electrical-transport measurements

All electrical measurements were conducted in a PPMS DynaCool cryostat (Quantum Design) with a base temperature of 1.8 K and a magnetic field up to 9 T. The conductance–temperature results were measured in a two-terminal configuration using a Model 7270 DSP lock-in amplifier to source a 1 mV AC voltage and to measure the current using a 13.333 Hz reference frequency, with a warm-up recipe at 1 K min^−1 being used. The other measurements were performed in a two-terminal configuration using a Keithley 2450 SourceMeter to source a DC voltage and measure current. The DC results were verified by AC measurements showing no substantial difference. The tunnelling junction dominates the final magnetoconductance owing to the high conductivity of the thin graphite contacts. The thickness of the prepared Ti/Au electrodes does not influence the experimental results. Note that a two-terminal configuration is widely used in 2D MTJs^26,27. The Shubnikov–de Haas oscillations of graphite electrodes observed in CrI₃ MTJs induced by out-plane field sweeping²⁵ were well suppressed in our experiments because we were only concerned with in-plane magneto-transport studies. The TMR ratio is given by (G_h − G_l)/G_l, in which G_l and G_h are the low conductance and high conductance, respectively. All field-direction-dependent transport measurements used the pre-calibrated local coordinates of the PPMS system and the angle is defined as ϕ. We used a unified method to mount the devices onto the electrical-transport puck by taking the ‘T’-distributed Ti/Au electrodes as the reference object. The ‘|’ electrode and the other two ‘—’ electrodes were aligned along the side edge and bottom edge of the puck by hand, respectively, resulting in the direction of the ‘|’ electrode being almost along the direction of the magnetic field when ϕ = 0°. This unified custom method helped us to examine the twist angles and locate the crystal axes of CrSBr flakes. We estimate that the uncertainty in angle is about ±5°.

DFT calculations for tunnelling mechanism

The electronic structure and the optimal interlayer separation are carried out based on DFT using the plane-wave projector augmented wave method⁴⁷, as implemented in the Vienna Ab initio Simulation Package (VASP)^48,49. We used the Perdew–Burke–Ernzerhof exchange-correlation functional⁵⁰ within the generalized gradient approximation (GGA). The electron–electron correlation effects beyond the GGA are taken into account using the Hubbard U correction through the GGA + U method⁵¹ with U value of 4.0 eV and Hund’s coupling energy J of 0.8 eV on the Cr 3d orbitals. For the self-consistent calculations, a plane-wave basis set with a plane-wave cutoff of 500 eV and a k-point mesh of 12 × 8 × 4 and 12 × 8 × 1 is used for the bulk CrSBr and thin films, respectively. The structural optimization is carried out using VASP maintaining the symmetry of the heterostructure. The positions of the atoms are relaxed towards equilibrium until the Hellmann–Feynman forces become less than 0.01 eV Å^−1.

Calculations of the complex band structure are performed using the non-equilibrium Green’s function formalism (DFT + NEGF approach)^52,53, as implemented in Synopsys QuantumATK⁵⁴. In QuantumATK, the nonrelativistic Fritz Haber Institute pseudopotentials are used with a single-zeta-polarized basis and the cutoff energy is set to 150 Ry. The spin-polarized GGA + U method is used in our calculations with the same U and J values on the Cr 3d orbitals as in VASP. A k-point mesh of 18 × 16 × 10 is used for bulk CrSBr. Periodic boundary conductions are assumed for the transverse direction (x–y plane) and open boundary conditions along the transport z direction. Experimentally measured lattice constants of a = 3.540 Å, b = 4.755 Å and c = 8.394 Å for the bulk CrSBr are assumed in the calculations³³. Figure 4d shows the distribution of the decay rates calculated at the energy of 0.1 eV below the conduction band minimum (CBM), that is, E_F = E_CBM − 0.1 eV, considering that CrSBr behaves as an n-type semiconductor. We also calculated the distribution of the decay rates in the 2D Brillouin zone at other energies and obtained TMR ratios ranging from about 200% to about 20,000% (Extended Data Fig. 8), which corroborates that our system inherently has a large TMR effect.

To take into account the effect of twisting on TMR, we assume that the in-plane wavevector \({\vec{k}}_{\parallel }\) characterizing the Bloch states of the two CrSBr layers, one twisted with respect to the other, is conserved in the tunnelling process, that is, coherent tunnelling. This approximation neglects the difference between the momentum and quasi-momentum and is applicable to a sufficiently thick vacuum barrier. When the top CrSBr monolayer is twisted with respect to the bottom monolayer, the wavevectors of the bottom monolayer remain unchanged as the reference, whereas the wavevectors of the twisted top monolayer are transformed as follows

Then the total transmission for the P and AP states as a function of twist angle θ can be obtained from

Equations (6) and (7) reflect the effect of twist on the transmission owing to the rotation of the in-plane wavevector but with collinear magnetizations. In the experiment, however, the magnetizations of the two CrSBr monolayers are non-collinear owing to the strong magnetocrystalline anisotropy. Therefore, the transmissions need to be calculated when the angle between the magnetizations of the two monolayers is θ (the same angle as for the structural rotation) and π − θ (magnetization of the top monolayer switched by 180°). Using the same idea as the Slonczewski model¹⁰, we quantize the spin of the top monolayer with respect to the spin axis of the bottom layer, obtaining

We can see that, at θ = 0°, equations (8) and (9) become equations (1) and (2), respectively. The TMR ratio for the twisted case as a function of θ is then given by

that is, equation (3). From the calculated distribution of the lowest decay rates as a function of \({\vec{k}}_{\parallel }=({k}_{x},{k}_{y})\) at E_F = E_CBM − 0.1 eV, we compute the distribution of \(\overrightarrow{{k}_{\parallel }^{{\prime} }}=({k}_{x}^{{\prime} },{k}_{y}^{{\prime} })\) as a function of different twist angle θ. Using equations (6), (7) and (3), we then calculate TMR as a function of twist angle θ.

Apart from non-collinear magnetizations, our model also considers coherent tunnelling, that is, conserved in-plane wavevector as discussed above. By contrast, in the Slonczewski model, \(G=\frac{{G}_{{\rm{P}}}+{G}_{{\rm{AP}}}}{2}+\frac{{G}_{{\rm{P}}}-{G}_{{\rm{AP}}}}{2}\cos \theta \), only the former is considered. Applying the experimental TMR ratio of the untwisted bilayer \(\frac{{G}_{{\rm{P}}}-{G}_{{\rm{AP}}}}{{G}_{{\rm{AP}}}}=\mathrm{1,050} \% \) (Extended Data Fig. 1d) in \(G=\frac{{G}_{{\rm{P}}}+{G}_{{\rm{AP}}}}{2}+\frac{{G}_{{\rm{P}}}-{G}_{{\rm{AP}}}}{2}\cos \theta \), we obtain the green curve in Fig. 4e.

To investigate the optimal interlayer separation, that is, the vdW gap, for the untwisted CrSBr bilayer, initially, in the AF configuration, we used the bulk lattice parameters and the bulk vdW gap of 2.78 Å at room temperature. By varying the vdW gap, we observed an optimal separation of around 2.39 Å at 0 K, which is 0.39 Å smaller than the bulk gap at room temperature, which is reasonable. Furthermore, we also explored the FM configuration, although it is not the experimental ground state of CrSBr. Notably, we found that the optimum vdW gap remained unchanged, indicating consistency across different magnetic configurations in CrSBr.

We further extended similar calculations to a twisted CrSBr structure. To facilitate the calculation, we used the following formula⁵⁵

and selected θ = 37.12°, a commensurate twist angle suitable for an orthorhombic structure determined by setting p = 1, q = 2 and 2δ as the angle indicated by the blue lines in Supplementary Fig. 7d. When analysing the interlayer separation versus energy for the twisted structure, we find that the vdW gap at 0 K is greatly increased to 3.36 Å. This outcome indicates that the twist operation induces a decoupling effect on the layers, thereby reducing the interlayer AF coupling and electron hopping. To further illustrate the weakening of interlayer coupling, we conducted projected density of states (DOS) calculations for a layer within this twisted CrSBr bilayer while maintaining the other layer unchanged (Supplementary Fig. 8e). For simplicity, we ensured that the magnetic atoms within each layer had the same spin direction. The projected DOS of an individual monolayer of CrSBr was also calculated with a similar spin configuration (Supplementary Fig. 8d).

DFT calculations for exchange interactions

The calculations were performed using a self-consistent Green function code Hutsepot based on the multiple-scattering theory, specially designed for semi-infinite systems real-space clusters^56,57. The calculations were performed within the GGA⁵⁰. Strongly localized Cr 3d electrons were treated using a combination of a self-interaction correction and the Slater transition-state methods as implemented within the multiple-scattering method^58,59,60,61. Similar results were also obtained using a GGA + U method^51,62. We used the maximum angular momentum cutoff of L_max = 3. The Brillouin zone integration was done using a tetrahedron method adapted for 2D geometry⁶³. The other parameters needed are similar to the previous calculation for the tunnelling mechanism. There were estimated exchange parameters J_ij entering the Heisenberg model

Here the unit vectors are e_i = S_i/|S_i| at site i, in which S_i is the localized spin moment. The exchange parameters were calculated using the magnetic force theorem as it is implemented within the multiple-scattering theory^64,65

in which \({G}_{\uparrow (\downarrow )}^{ij}\) is a Kohn–Sham Green’s function between the sites i and j for the spin-up (spin-down) channel and Δ_i is a magnetic interaction vertex function. Supplementary Fig. 9 shows the schematic of the magnetic exchange interactions. J₁ to J₇ denote the intralayer exchange interactions and J_z1 and J_z2 denote the interlayer exchange interactions. The calculated exchange parameters are shown in Extended Data Table 2. The first column shows that the magnetic interaction within each monolayer of CrSBr bulk is strongly ferromagnetic (J₁ to J₇ are mainly positive) and, between the neighbouring layers, it is antiferromagnetic (J_z1 < 0, J_z2 is small). The magnetic coupling is mainly mediated by superexchange: within a monolayer through S sp orbitals (S has an induced magnetic moment of 0.08μ_B) and between two monolayers through Br sp orbitals (with the induced magnetic moment of 0.03μ_B). Cr and Br atoms form a long linear bond under the angle of 180°, Cr–Br–Br–Cr (Supplementary Fig. 9b), which is responsible for an AF interaction (J_z1 = −0.25). In the same time, J_z2 remains relatively small in magnitude because: (1) the distance between Cr moments over the vdW gap is large for a direct coupling and (2) the bonds between Cr and the neighbouring anions (S and Br) are rotated by a large angle (the bond strength between two monolayers in this direction is rather weak; see Supplementary Fig. 9b). The theoretically obtained Néel temperature 132 K using a random-phase approximation was found to be in good agreement with the experiment. We also calculated the exchange parameters for CrSBr monolayer and bilayer as shown in Extended Data Table 2. The values near the exchange parameters of the CrSBr bulk prove that FM coupling is robust to the monolayer and AF coupling is robust to the bilayer.

Then we investigated the exchange interactions in twisted CrSBr bilayers with twist angles of 45° and 90°, respectively. To calculate exchange parameters in a twisted bilayer, the Green function in equation (12) was calculated using the transformation

in which \({G}_{{{LL}}^{{\prime} }}^{{ij}}(E)\) are Green function matrix elements and \({U}_{{LL}^{{\prime} }}({s}_{i}\,;E)\) is the transformation matrix for a displacement s_i (refs. ^66,67). The calculated exchange parameters are also shown in Extended Data Table 2. The intralayer coupling is changed only marginally, whereas the interlayer coupling was substantially suppressed. The reason for this is that the Cr–Br–Br–Cr bonds are weakened on twist operation: Cr and Br are not lying on one line but are tilted; see Supplementary Fig. 9c,d. This leads to a marked decrease of superexchange between Cr atoms through Br sp orbitals across the vdW gap.

Overall, the DFT results show a negligible interlayer magnetic coupling at the twisted interface and a strong intralayer FM coupling, which support that both quasi-parallel and quasi-antiparallel spin alignments are stable at ZF, as shown in Fig. 1f. Therefore, such bistable spin alignments at the twisted interface can be operated by sweeping external field, being well described by Stoner–Wohlfarth model calculations (Supplementary Note 3). Moreover, the DFT results rather establish a strong interlayer AF coupling at the untwisted interfaces, which means that the net magnetization is zero in the twisted bilayer/bilayer MTJs.

Strong field behaviours of twisted bilayer/bilayer MTJs

For the 35° twisted bilayer/bilayer MTJ, forward and backward sweeping between −1.5 T and 1.5 T is applied in the a–b plane (Extended Data Fig. 2a,b). Compared with Fig. 2e and Extended Data Fig. 1a, a twofold rotational symmetry is preserved, but the smooth relationship between the orientation of the field and the saturation field in the case of single bilayer MTJ is broken, as shown in Extended Data Fig. 2a,b. A distinct feature is that two tiny humps appear in the vicinity of the two a axes. The angle between the positions of the two humps precisely matches the twist angle. These results demonstrate that the original magnetic anisotropy in the intrinsic CrSBr bilayer is robust, irrespective of the twisted alignment, and strongly support that each bilayer in the twisted structure is largely independent and magnetically decoupled from each other. Apart from the inherent physical properties of CrSBr, the decoupling may be a result of the large twist angle as distinct from small-angle-twisted moiré superlattices (usually less than 10°)^{46,68,69,70,71,72}, in which the interlayer coupling mechanism is prevalent. To further experimentally confirm whether there is substantial interlayer magnetic coupling in our system, we also performed control experiments on a 30° twisted CrSBr bilayer/hBN monolayer/CrSBr bilayer MTJ (Supplementary Fig. 10), together with the twisted monolayer/hBN monolayer/monolayer MTJ (Extended Data Fig. 7). The inserted hBN layer would reduce or eliminate the hypothetical coupling at the twisted interface, but we still observe ZF NV in both cases, which affirms the aforementioned independence and interlayer decoupling in line with the DFT calculations. Further experimental evidence is that the symmetry of the tunnelling current curve is mirrored, as well as the values of the two ZF tunnel currents being maintained after a large field application (Extended Data Fig. 2g,h) owing to the spin configuration being flipped to its time-reversal copy in the pinned bilayer, which is reproduced by Stoner–Wohlfarth model calculations in Supplementary Note 3. This is also the reason why the ZF NV is unstable using strong field sweeping (Extended Data Fig. 2c–f).

Distinct temperature-dependent behaviours between twisted and untwisted interfaces

In the twisted bilayer/bilayer MTJ (Extended Data Fig. 4k), the high-temperature measurements demonstrate that the ZF-TMR calculated from the conductances of ZF-on versus ZF-off resulting from the twisted interface is robust up to close to T_N. By contrast, the TMR ratio calculated from the conductance at 9 T and ZF-on quickly decays with increasing temperature. Note that the twisted bilayer/bilayer has three interfaces, that is, one twisted interface and two untwisted interfaces (within each bilayer). Because the spin alignments are parallel at all three interfaces of the twisted bilayer/bilayer MTJ at 9 T and ZF-on corresponds to quasi-parallel spin alignment at only the twisted interface (the right side of Fig. 1f), the TMR from 9 T versus ZF-on mainly originates from the two untwisted interfaces. Thus, these distinct temperature-dependent behaviours strongly imply the uniqueness of the twisted interface in contrast to the untwisted interface. We then confirmed this discrepancy by comparing an untwisted bilayer MTJ and a twisted bilayer MTJ (Extended Data Fig. 6f). Note that the marked differences in the temperature-dependent variations of TMR for the twisted and untwisted interfaces are independent of the applied bias (Supplementary Figs. 11 and 12).

The reason for this is ascribed to the decoupling mechanism at the twisted interface. In our 2D MTJs, the whole CrSBr stack plays the role of a tunnelling barrier. As a semiconductor, CrSBr is insulating at low temperatures, which well serves as the tunnelling barrier, whereas CrSBr becomes conductive with increasing temperature (Fig. 2c), resulting in degradation of the tunnelling barrier, which can thereby naturally explain the fast decay of the TMR in the untwisted MTJs. Nevertheless, this cannot account for the slow decay behaviour in the twisted MTJs. We have already noted the important difference between the untwisted and twisted interfaces from DFT calculations of the interlayer exchange interactions (Extended Data Table 2). Now we consider the energy-band structure. Supplementary Fig. 8a,b shows the energy-band structure of the untwisted CrSBr bilayer, in which the energy band of the top and bottom monolayers (Supplementary Fig. 8c) hybridizes or overlaps well, which can be understood from a more intuitive picture. Supplementary Fig. 8d shows the projected DOS calculated from an individual CrSBr monolayer, from which it is found that Cr d orbitals dominate the DOS near the CBM. For the case of two untwisted monolayers, the Cr d orbitals in each monolayer are well aligned to allow for an effective interlayer orbital interaction by means of the intervening Br atoms, resulting in superexchange across Cr–Br–Br–Cr. Thus, increasing temperature can effectively enhance the electron interlayer hopping across the vdW gap. Such thermal activation is spin-independent and we surmise that this mechanism progressively dominates the vertical conductance of the untwisted bilayer on warming. However, because the Cr d orbitals are highly spatially anisotropic, rotating one monolayer with respect to the other should substantially reduce the effective interlayer orbital interaction⁷³. Supplementary Fig. 8e shows the projected DOS of the top monolayer of a twisted bilayer. Comparing Supplementary Fig. 8e with Supplementary Fig. 8d, the slight difference indicates that interlayer orbital interaction with the bottom monolayer is indeed suppressed at the twisted interface in accordance with the negligible interlayer magnetic coupling suggested by DFT calculations. Accordingly, we suggest that the thermally activated spin-independent interlayer hopping is less in the twisted MTJs and the spin-dependent tunnelling mechanism still contributes part of the vertical conductance. Experimental evidence for this is a slower increase of conductance with warming (Extended Data Fig. 6b) compared with the conductance–temperature curve of the untwisted MTJ (Fig. 2c). On the other hand, we find that the vdW gap becomes physically thicker in the twisted MTJs (Supplementary Fig. 13). DFT calculations confirm that a twisted structure with a thicker vdW gap is energetically favourable (Supplementary Fig. 7). Therefore, the thicker vdW gap can still serve as a robust tunnelling barrier (even though the CrSBr sheets themselves become more conducting) up to close to T_N, giving a robust weakly temperature-dependent TMR in the twisted MTJs.