Tunneling current-controlled spin states in few-layer van der Waals magnets

Effective control of magnetic phases in two-dimensional magnets would constitute crucial progress in spintronics, holding great potential for future computing technologies. Here, we report a new approach of leveraging tunneling current as a tool for controlling spin states in CrI3. We reveal that a tunneling current can deterministically switch between spin-parallel and spin-antiparallel states in few-layer CrI3, depending on the polarity and amplitude of the current. We propose a mechanism involving nonequilibrium spin accumulation in the graphene electrodes in contact with the CrI3 layers. We further demonstrate tunneling current-tunable stochastic switching between multiple spin states of the CrI3 tunnel devices, which goes beyond conventional bi-stable stochastic magnetic tunnel junctions and has not been documented in two-dimensional magnets. Our findings not only address the existing knowledge gap concerning the influence of tunneling currents in controlling the magnetism in two-dimensional magnets, but also unlock possibilities for energy-efficient probabilistic and neuromorphic computing.


Supplemntary Note 1: Tunneling magnetoresistance in the linear response regime
To calculate the tunneling magnetoresistance for our model it is sufficient to consider the two insulator layers only.In this case we have (Σ  < )  () = 2Γ (  1 (ℏ)  2 (ℏ) ) where  1,2 ± = ℏ −    + Δ +  1,2  +   ± Γ.It is then straightforward to get ). (2) In the present case the (electron number) current operator is The tunneling current for each spin and orbital is therefore where ′ ≡ ℏ,   () is the density of states in each insulator layer, and in the 2nd line we have omitted the  subscripts.Moreover, (5) where we assumed  = 0 and  is much smaller than the band gap.When  1 =  > 0 and  2 = 0, Eq. 4 becomes where we have further assumed that the density of states   () ≈   when  ∈ [−  2 ,  2 ],  being the band width, and 0 otherwise.Finally, if both  and  are much smaller than the band gap, the integrand is approximately a constant.We thus obtain which means the tunneling conductance   ≡   / is In the case of CrI3 since the gap is between the two lowest bands having the same spin, the chemical potential   = 1 2 (− − Δ −  + Δ) = −.We can thus obtain, when the magnetizations of the two layers are parallel: while when they are antiparallel One can further obtain the tunneling magnetoresistance: For moderate values of Δ/, e.g.Δ/ = 0.7,0.8,0.9, we find that TMR = 192%, 109%, 60%.We therefore chose Δ/ = 0.8 to get the results in the main text so that the TMR is comparable to previous experimental results.

Supplementary Note 2: Vanishing spin current in the antiferromagnetic state of the model without asymmetry
In this section we show that the tunneling spin current exactly vanishes in the antiferromagnetic state of our model when the two insulator layers are identical except for  1 = − 2 = .This is most easily seen by using the standard Landauer formula applied to our model: where   , are given in Methods, and below we will omit the subscript  for brevity.The tunneling currents obtained from this formula are consistent with that from the average inter-layer current discussed in Methods.This is a result of the fact that our model does not have other relaxation mechanisms except that through the lead self-energies and therefore the tunneling is fully coherent.When there is no asymmetry between the two insulator layers, in the antiferromagnetic state we have (  , ) ↑ being equivalent to (  , ) ↓ upon inverting the whole system: (  , ) ↑ =  † (  , ) ↓ ,   = ( 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 Therefore ( ↑ , )  = ( †  ↓ , )  = ( ↓ , ) ̅ ̅ = ( ↓ , ) ̅ ̅ (14) where in the last step we have used   =    .Eq. 14 indicates that, for example, ( ↑ , ) 14 = ( ↓ , ) 14 , ( ↑ , ) 11 = ( ↓ , ) 44 .Using this result we can obtain, along with Eq. 12: Namely, the tunneling spin current exactly vanishes in the antiferromagnetic state.This result also suggests that as long as the time-reversal plus spatial inversion symmetry of the four-layer model in the antiferromagnetic state is broken, the tunneling spin current should generally be nonzero.

Supplementary
Here   is the spin relaxation time.Below we first show that the lead self-energies result in an effective   and hence the spin accumulation in the present model.
The spin current for the present model is defined as For convenience we use the Landauer formula Eq. 12 for   .Therefore We next directly calculate the nonequilibrium spin density.The spin densities in the top and bottom layers are where the relevant matrix elements of   < can be further expressed as (below we drop the subscript  for brevity) However, the spin densities in the above expressions also include equilibrium contributions, since unlike the current, 〈  〉 does not vanish in equilibrium.Moreover, the local spin and charge densities can also change due to shifting the chemical potentials in all leads by the same amount, which is clearly an equilibrium effect.Such an equilibrium effect is also existent when  1 ≠  2 and is intertwined with the tunneling current induced nonequilibrium spin accumulation.Although generally speaking in a tunneling setup it may not always be possible to separate equilibrium and nonequilibrium contributions unambiguously for quantities that are nonzero in equilibrium, for the present model where there is only one metal layer attached to each lead it is sensible to define the nonequilibrium contributions as This is equivalent to changing Eq.21 into As a result Therefore the leads self-energies correspond to a spin-relaxation time Because spin-relaxation in our model occures in the leads,   model here is equal to the escape time from the device and therefore related to the contact resistance.In a more realistic model which incorporates spinrelaxation within the device, these two quantities can vary independently.
Based on the above discussion, we calculate the nonequilibrium spin accumulation in the top and bottom metal layers in this work by first using Eq.24 and then rescaling them by   exp /  model , with   exp (or simply   used in the main text) the experimentally determined spin relaxation time in graphene.
2. The stochastic and unidirectional switching events occur ideally in different ranges of current densities.At high current densities significant heating leads to deviation of the dynamic magnetic states in the bilayer from the static SAP and SP states.In such a case conventional spin-transfer torque is effective and indeed the SP state always seems to be favored over the SAP state at the high-current end of the stochastic switching regions (see e.g., Fig. 3a in the main text).However, different domains in general have different high-and low-bias ranges thus defined, and the "anomalous" SAP to SP switching at negative bias can originate from the low-bias end of the stochastic switching range of a different domain.In addition, we also found a similar SAP to SP switching feature in the negative-bias range between 0.58 T and 0.575 T in Fig. S6a, which as the magnetic field keeps decreasing merges into the stochastic switching region at around 0.575 T.
the working devices as the graphene/CrI3/graphene tunneling junctions with nonlinear I-V characteristics (characteristics of the tunneling behavior) and clear magnetic field dependent tunneling magnetoresistance.
Supplementary Note 8: Estimate the switching current density of the unidirectional spin-state transition We calculated the switching current density by dividing the tunneling current by the estimated area of a single domain in CrI3 by directly analyzing the transport data in Fig. 2a.To extract the area of the magnetic domain and the current responsible for the spin state transition, we assume the Ohm's law is valid, and the resistances of the magnetic domains with SP and SAP states connected in parallel in the circuit as shown in Fig. S5. Figure S5a-c depicts the simplified circuits of the 2L CrI3 tunnel junction device displaying a layer SP state, a layer SAP state, and two magnetic domains with both SP and SAP states, respectively.As demonstrated in Fig. 2b of the main text, near the SAP to SP transition region, the layer magnetic state of the 2L CrI3 is SAP when the magnetic field is below 0.2 T and shifts to SP when it surpasses 0.7 T.
Furthermore, Figure 2b reveals that the current-induced transition from the SAP state (marked as ①) to the SP state (marked as ②) can be interpreted as the SAP state (with two magnetic domains in the layer) flipping to SP (with only one magnetic domain now).It's worth mentioning that there may be one or two additional small domains.However, their influence on the resistance change is minimal compared to the SAP to SP transition highlighted in the inset of Fig. 2a of the main text.Now, we can use the data in Fig. S5d to estimate the area of the magnetic domain with the SAP state.When the magnetic field is above 0.7 T, the layer magnetic state is SP, and we have   =     , where RSP is the resistance of the junction, A is the junction area (~ 0.5 µm 2 ), L is the thickness of the sample, and   is the corresponding resistivity.At 0.7 T (see Fig. S5d), we have    =    = 0.5490  × 0.5  2 = 0.2745  •  2 .Similarly, at 0.2 T (Fig. S5d), the layer magnetic state is SAP, and we have   =     , where RSAP and   are the resistance and resistivity of the tunnel junction, respectively.Then,    =     = 0.603  × 0.5  2 = 0.3015  •  2 .In the case of the intermediate state with one SP domain and one SAP domain (Fig. S5c) in the 2L CrI3, the total resistance R' of the intermediate state is , where ′  and ′  are the resistances of the SP and SAP domains, respectively.Then, we have , where   is the area of the domain with an SAP state.From Fig. S5d, we know  ′ = 0.573 MΩ, resulting in the area (  ) of the magnetic domain with an SAP state of ~0.235 µm .From Fig. S5d, we know at state ①,  ′ = 0.573 MΩ, resulting in the area (  ) of the magnetic domain with an SAP state is ~0.235 µm 2 .At state ② and  ′ = 0.556 MΩ, the estimated area of the magnetic domain is ~ 0.072 µm 2 .Thus, the area of the magnetic domain that got switched from ① to ② is ~ 0.163 µm 2 , which is comparable with the areas of the magnetic domains of CrI3 thin layers captured by single-spin microscopy, as referenced in Ref. 35.Next, we calculate the switching current,   , which is responsible for the SAP to SP transition at the bias current I of 3 µA.From Fig. S5d, we have   = ′  ′  +′  = 1.18 µA.Thus, the switching current density (

𝐼 𝑆𝐴𝑃 𝐴 𝑆𝐴𝑃
) is ~724 A/cm 2 , being around three orders of magnitude lower than the values reported in previous studies employing SOT.We note that the transition between SAP and SP states can be realized at a much lower current of ~100 nA, and the switching current density could be even lower.

Fig. S1|
Fig. S1| Cooling history-dependent magnetic domain formation in 2L and 4L CrI3.Tunneling resistance as a function of an applied magnetic field with different cooling cycles for a-c a 2L CrI3 tunnel junction device and d-f a 4L CrI3 tunnel junction device.All measurements were performed at T = 1.5 K.One cooling cycle represents warming up the device to room temperature and cooling it down to the base temperature again.

Fig. S2|
Fig. S2| Temperature dependence of tunneling magnetoresistance for 2L and 4L CrI3.Tunneling resistance as a function of an applied out-of-plane magnetic field at different temperatures at a constant bias of 350 mV for a a 2L CrI3 and 300 mV for b a 4L CrI3, respectively.

Fig. S3|
Fig. S3| Bias voltage-dependent coercivity of a 2L CrI3 (Device 2). a Tunneling current as a function of both bias voltage and applied magnetic field of a 2L CrI3, where the magnetic field was swept from -1 to 1 T. b Selected line profile from the white-dashed line indicated in (a).The measurements were performed at T = 1.5 K.

Fig. S4|
Fig. S4| Bias voltage-dependent coercivity of a 4L CrI3. a and b Tunneling current as a function of both bias voltage and magnetic field of the 4L CrI3, where B was swept from 1.9 to 1.55 T (a) and from 1.55 to 1.9 T (b).c Selected line profiles from the white-dashed line indicated in (a) and (b).The measurements were performed at T = 1.5 K.

Fig
Fig. S5| Simplified circuits for estimating the size for magnetic domain.The circuits for 2L CrI3 tunnel device with a layer SP state, b layer SAP state, and c two magnetic domains with both SP and SAP states.dTunneling resistance as a function of the magnetic field at a bias current of 1 μA (top panel of Fig.2a).

Fig. S6|
Fig. S6| Magnetic field dependent I-V characteristics for a 2L CrI3 at T = 1.5 K. a Voltage as a function of applied current at different magnetic fields ranging from 0.5 to 0.6 T with a step size of 2 mT.b Representative I-V curves measured at magnetic fields of 0.588, 0.534, and 0.510 T, respectively.Insets are zoom-in curves from the corresponding dashed green boxes.

Fig. S7|
Fig. S7| Magnetic field dependent I-V characteristics and current coercivity for a 2L CrI3. a Voltage as a function of applied current at different magnetic fields ranging from 0.5 to 0.6 T with a step size of 2 mT.b Representative I-V curves measured at magnetic fields of 0.514T.Inset is zoom-in curves from the dashed black box.c B-field-dependence of the current coercivity for the hysteresis loop around 3000 nA, as shown in the inset of (b).The Ilow and Ihigh are the current of the red and blue curve at the transition points, respectively.∆ = Ihigh -Ilow.All the data were measured at 1.5 K.

Fig. S8|
Fig. S8| Time snapshot of voltage for an applied DC current in the fluctuation region with a time scale of 5 s of a 2L CrI3.The measurement was performed at B = 0.550 T and T = 1.5 K.

Fig. S9|
Fig. S9| Magnetic field dependent I-V characteristics for a 2L CrI3 at T = 10 K. a Voltage as a function of applied current at different magnetic fields ranging from 0.5 to 0.6 T with a step size of 2 mT.b Representative I-V curves measured at magnetic fields of 0.588, 0.564, and 0.518 T, respectively.Insets are zoom-in curves from the corresponding dashed green boxes.

Fig. S10|
Fig. S10| Magnetic field dependent I-V characteristics for a 2L CrI3 at T = 20 K. a Voltage as a function of applied current at different magnetic fields ranging from 0.5 to 0.6 T with a step size of 2 mT.b Representative I-V curves measured at magnetic fields of 0.588, 0.550, and 0.500 T, respectively.Insets are zoom-in curves from the corresponding dashed green boxes.

Fig. S11|
Fig. S11| Magnetic field dependent I-V characteristics for a 2L CrI3 at T = 30 K. a Voltage as a function of applied current at different magnetic fields ranging from 0.5 to 0.6 T with a step size of 2 mT.b Representative I-V curves measured at magnetic fields of 0.586, 0.572, and 0.546 T, respectively.Insets are zoom-in curves from the corresponding dashed green boxes.

Fig.
Fig. S12| Current-driven magnetization reversal in a 4L CrI3. a Voltage as a function of applied current of a 4L CrI3 measured at B = 1.720 T. b Tunneling resistance as a function of an applied out-of-plane magnetic field at different bias currents.The initial magnetic states are prepared at B = 1.720T, then the current is ramped to 1 μA and then back to the target currents as indicated in each graph.The circled numbers in (a) and (b) indicate the corresponding initial spin states.All the measurements were performed at T = 1.5 K.

Fig
Fig. S13| I-V characteristics of a 4L CrI3 tunneling device.a Optical image of a 4L CrI3 tunneling device.b I-V characteristics of the 4L CrI3 tunneling device measured at B = 0 and 2 T. The measurements were performed at T = 1.5 K.

Fig. S14|
Fig. S14| Bias voltage-dependent coercivity in a 4L CrI3. a and b Tunneling current as a function of magnetic field of the 4L CrI3 tunneling device measured at both positive (a) and negative (b) bias voltages.All the measurements were carried out at T = 1.5 K.

Fig. S15|
Fig. S15| Magnetic field dependent I-V characteristics for a 4L CrI3.a-l Voltage as a function of applied current at different magnetic fields near the layer-magnetic phase transition region.All the measurements were performed at T = 1.5 K.

Fig. S16|
Fig. S16| Stochastic switching of a four-layer (4L) CrI3. a Voltage as a function of applied current measured at B = 1.715 T. The inset zoomed from the dashed black square shows the current-driven spin state transition.b V versus I curve zoomed from the dashed green square in (a), showing voltage fluctuations between three spin states.c Left panels show time snapshots of voltage for constant currents of -760, -730, -700, -530 and -520 nA (from bottom to top).Right panels show the corresponding histograms of voltage distribution with a sampling time of 600 s.All the measurements were performed at T = 1.5 K.

Fig
Fig. S17| Signatures of unidirectional magnetization reversal and stochastic switching in a 5L CrI3 magnetic tunneling device.a I-V curves were taken at T = 1.5 K, B = 0.61 T and b at T = 1.5 K, B = 0.586 T. Black (red) line corresponds to forward (backward) sweeping.

Fig. S18| 1 𝛤
Fig. S18| Tunneling electron and spin currents (a, c, e) and nonequilibrium spin accumulation (b, d, f) in the SAP with different types of asymmetries.a, b The top insulator layer ( 1 < 0) has a smaller orbital splitting (Δ) than that of the bottom layer by 0.1 eV.c, d The bottom insulator layer ( 2 > 0) has a higher chemical potential (  ) than the top layer by 0.1 eV, mimicking the effect of a perpendicular electric field.e, f The hopping between the bottom insulator layer ( 2 > 0) and the bottom metal layer is larger than that between the two top layers by 0.01 eV.