Spin‐Torque Memristors Based on Perpendicular Magnetic Tunnel Junctions for Neuromorphic Computing

Abstract Spin‐torque memristors are proposed in 2009, and can provide fast, low‐power, and infinite memristive behavior for neuromorphic computing and large‐density non‐volatile memory. However, the strict requirements of combining high magnetoresistance, stable domain wall pinning and current‐induced switching in a single device pose difficulties in physical implementation. Here, a nanoscale spin‐torque memristor based on a perpendicular‐anisotropy magnetic tunnel junction with a CoFeB/W/CoFeB composite free layer structure is experimentally demonstrated. Its tunneling magnetoresistance is higher than 200%, and memristive behavior can be realized by spin‐transfer torque switching. Memristive states are retained by strong domain wall pinning effects in the free layer. Experiments and simulations suggest that nanoscale vertical chiral spin textures can form around clusters of W atoms under the combined effect of opposite Dzyaloshinskii–Moriya interactions and the Ruderman–Kittel–Kasuya–Yosida interaction between the two CoFeB free layers. Energy fluctuation caused by these textures may be the main reason for the strong pinning effect. With the experimentally demonstrated memristive behavior and spike‐timing‐dependent plasticity, a spiking neural network to perform handwritten pattern recognition in an unsupervised manner is simulated. Due to advantages such as long endurance and high speed, the spin‐torque memristors are competitive in the future applications for neuromorphic computing.

The MTJ films were deposited by a Singulus TIMARIS 200 mm magnetron sputtering machine at a base pressure of 3.75×10-9 Torr.

S1.1 Saturation magnetization and perpendicular anisotropy energy
Both the in-plane and out-of-plane (perpendicular) hysteresis loops of the film used to fabricate the magnetic tunnel junction (MTJ) were measured via vibrating sample magnetometry (VSM), as shown in Figure S1a. Both the free layer (FL) structure and the reference layer exhibited good perpendicular magnetic anisotropy (PMA). Furthermore, the dual FLs were globally well ferromagnetically coupled.
After producing a film (FL-film) with the same structure as the FL structure of the MTJ stack, we characterized the in-plane ( Figure S1b  The perpendicular anisotropy constant KU of the FLs was estimated based on the in-plane field hysteresis loop: the critical field that rotates the magnetization of the FLs to the in-plane direction is considered as the effective anisotropy field 0 K,eff . The uniaxial anisotropy energy was calculated as follows: [1] = 1 2 0 , + 1 2 0 2 . At room temperature, the effective anisotropy field 0 K,eff of the FL was measured to be 250 mT, and the anisotropy KU was estimated to be 7.5×105 J/m3.

S1.2 Characterization of the FL-film via NV center in diamond
Stray fields above the FL-film is scanned using NV center in diamond to check the homogeneity of the magnetization of the film. The NV center is created by implantation of 30 keV 2 + ions into [100] bulk diamond and annealing for 2 hours at 1000 ℃ in vacuum.
Therefore, the depth of the NV center is estimated to be 20.7 nm.
As a comparison, we numerically calculate the fluctuation of the stray field above a sample caused by an eventual defect, where the magnetization is locally changed. We suppose a magnetic film with homogenous magnetization = 1.0 × 10 6 A/m , except for one 10 nm×10 nm defect in the center where the magnetization reduces by 10%, as shown in Fig S2a. The thickness of the film is 1.8 nm and the magnetization is saturated in the z-direction. The NV center is placed 20 nm above the sample plane and the NV quantization axis is in the X-Z plane, tilting 35° from the sample plane. Micromagnetic simulation software MuMax3 [2] is used to calculate the stray field projected along the NV quantization axis in the detection plane. The area considered in the simulation is 150 nm×150 nm, with 20×20 periodic repetitions to minimize the effect of boundary conditions. The simulation result ( Figure S2b) shows that a fluctuation of stray field caused by such an isolated defect is in the order of 300 μT. Whereas, Figure S2 Numerical evaluation of the stray field detected by the NV centre. a) Schematic diagram of simulation configuration. b) Distribution of the stray field in the detection plane, projecting along the NV quantization axis.
the experimentally detected stray fields via the NV center are in the order of 80 μT. This comparison proves that the magnetization of the FL of the studied film is very homogeneous.
It should be noticed that the bias of stray field shown in Figure S2b is caused by the finite dimension of the simulated area, which is much smaller than the experimentally measured sample. when driven by a perpendicular field under the application of a varying in-plane field. [3,4] High perpendicular fields ⊥ are used so that the DW moves beyond the thermally activated creep regime. Pinning effects enhanced by the application of in-plane fields [5] have no influence in this regime and the measured DMI is believed to be more accurate. [4] The results are shown in Figure S3. In each measurement, the perpendicular remains unchanged while the in-plane field varies step by step. The in-plane field under which the DW velocity reaches a minimum value can be regarded as the effective DMI field HDM. Accordingly, the strength of the overall DMI can be estimated as = 0 √ ⁄ . Here, ex is the Heisenberg exchange stiffness and eff is the effective anisotropy constant.

S1.4 Ruderman-Kittel-Kasuya-Yosida (RKKY) exchange through the W spacer layer
S. Parkin measured the RKKY exchange of two magnetic layers separated by a W layer. [6] The structure used was Si/Cr (3.5 nm)/[Co (1.5 nm)/W (x nm)]16/Cr (2 nm), with in-plane magnetic anisotropy in the Co layer. The following results were reported: the first antiferromagnetic coupling peak occurs at tW=0.55 nm, and the antiferromagnetic coupling strength is 0.03 mJ/m2 at this peak. The thickness of W spacer that corresponds to this first antiferromagnetic region is approximately 0.3 nm.
We fitted these data using the well-known RKKY law: [7] ( ) = • (2 ) (S1) where A is a parameter related to the magnitude of the coupling strength and kF is a parameter related to the oscillation period. The fitting results are shown in Figure S4. From this fit, we can conclude that the ferromagnetic coupling between the two magnetic layers is Jex=0.6 mJ/m 2 at tW=0.2 nm and drops to zero at tW=0.43 nm.

Section 2. Calculation of the domain wall surface energy
In this calculation, we suppose that the two FLs are mirror-symmetric with respect to the spacer layer. The tilt angle of the magnetization in the center of a DW ( ⃗⃗ DW ) with respect to the transverse direction is θ (-θ) in the LwFL (UpFL), where ∈ [0, /2]. The total magnetic energy of the DW is [7,8] = 0 + (1 − cos 2 ) ∆⁄ − 2 sin + (S3) where 0 = 4 √ is the surface energy of a Bloch DW with θ=0, [9] as described in the first case in Figure 5b, but without considering the DMI energy, Aex is the exchange stiffness, Keff is the effective anisotropy field, ∆ is the DW width, M is the thickness of a single FL, and is the energy associated to the dipole-dipole interaction.
The local dipole-dipole interactions between the DWs in the two FLs promote the formation of a chiral vortex DW because this type of DW configuration has a lower magnetic static energy than a coupled Bloch DW. [10] Since the theoretical calculation of this energy is very difficult, here, we perform a numerical calculation using MuMax3. [2] The simulated structure is shown in Figure S5. in a film with PMA [10] ). The saturation magnetizations of the two magnetic layers were both 1.0×10 6 A/m, and the perpendicular anisotropy energy was 7.5×10 5 J/m 3 . The exchange stiffness in each layer was Aex=1.3×10 -11 J/m. A coupled Bloch DW configuration was created as the initial state. Then, the DW center magnetization ⃗⃗ was tuned by means of a locally applied external field. The angle of ⃗⃗ in the LwFL (UpFL) θ with respect to the y-axis was varied from 0 to π/2 (0 to -π/2). The total demagnetizing energy for the simulated system, , was saved during the simulation.
Then, the DW energy per unit area associated with the dipole-dipole interactions was calculated as follows: where 0 is the demagnetizing energy of the whole simulated system when θ=π/2, which is taken as the reference level of the demagnetizing energy; =2 nm is the total magnetic layer thickness; and w=16 nm is the width of the simulated structure, i.e., the length of the DW.
The variation in as a function of θ is plotted in Figure S6, which can be fitted by ≈ (1 + cos 2 ). For this simulated configuration, the DW energy associated with the dipole- dipole interactions decreases by more than 1.6 mJ per unit area when the DW structure changes from the coupled Bloch configuration to the chiral vortex configuration. This is an important factor for the formation of chiral vortices.
After substituting with (1 + cos 2 ) in Eqution S3, we get, where A≈0.78 mJ/m2 in our case. By minimizing the DW surface energy using we obtain the value of θ at equilibrium: For eff = 1.25 × 10 5 / 3 , as measured via vibrating sample magnetometry (VSM) for the FLs, ex = 1.3 × 10 −11 J/m, [11] Δ=5 nm, [8] tM=1 nm, and D=±0.5 mJ/m2, the variations of θ and the DW surface energy as functions of the interlayer coupling strength Jex is obtained, as plotted in Figure 5c in the main text.

Section 3. Effect of the demagnetizing field on the stability of an intermediate state
The global demagnetizing energy of an FL is reduced if the FL is partially switched because the net magnetization approaches zero. In other words, a demagnetizing field Bdemag helps to stabilize a DW in an FL disc. Here, we use the concept of the effective magnetization current to evaluate the influence of such a demagnetizing field on the DW dynamics in the FL. [12,13] For  In the beginning, a DW was set at the left of the film as the initial state. Then a 10 mT perpendicular field was applied permanently. We observed that the DW moved to the right and then entered the area where RKKY interactions lost. At the same time, the structure of the DW changes from Bloch type to chiral vortex. As calculated before, the chiral vortex wall can minimize both the DMI energy and demagnetizing energy. Therefore, this structure served as an energy well and robustly pinned the DW. No further DW motion occurred in the following.
The magnetic state of the film was saved during the simulation, and two frames of them can be found in Figure S8 and the video can be found in the Supplementary material Movie S1 (separate file). Movie S1 (separate file). Movie showing the structural transformation of a DW when it moves from a region with strong interlayer coupling to a region with weak interlayer coupling, and the DW pinning effect induced by the formation of the chiral vortex. Detailed simulation parameters are given in section S4 of supplementary material S1.

Section 5. Model for the memristive behavior
For both the P to AP and AP to P switching process in our device, the variation of the resistance with time is difficult to be accurately described with a theoretical model. However, it can be phenomenologically fitted by an exponential function, with a characteristic time τ, as shown in Figure 6a in the main text. ln is found to be linear to the magnitude of applied voltage, as shown in Figure S9. In detail, for the P to AP switching, the change of device resistance with time under a certain voltage V can be described as, where ln + = + + + . Here, we defined the voltage inducing P to AP switching as positive, and + represents the characteristic time for P to AP switching. Similarly, for AP to P switching, (S10) Furthermore, considering the stochastic variation of the resistance, which is caused by stochastic depinning of domain walls, a noise factor is multiplied on the δ .
where is a random number with a normal distribution generated for the ℎ operation. Note that if is positive, multiplicative noise is also positive and the resistance R changes with a finite amplitude; if is negative, equals zero and R does not change after the operation, corresponding to a domain wall pinning event. By carefully adjusting the mean value and standard deviation of , the stochasticity of resistance change caused by devices' variation or domain wall pinning events can be adequately mimicked in the simulation. Figure 6a & b in the main text give examples of the memristive switching process obtained using the above model for fixed voltage. It is to note that we found it difficult to obtain a complex switched state under relatively low voltage during the AP to P switching process in our devices. This phenomenon is considered in the simulation and the value of in Equation S10 used in the simulation is set to be 115 Ω.

Section 6. Spike neural network simulations
In this section, we will give a detailed description of SNN simulation based on the spintorque memristor device, which could be divided into the following four steps.
The first step is the pretreatment of the input images using the receptive field function to make the input data closer to what a retina neuron would perceive. [14] As can be seen from FigureS10, where the input is a 20×20-pixel black-and-white image, we constructed an oncentered 5×5 receptive field array in which stimulation leads to response of a particular sensory Figure S9. The relationship between the characteristic time constant τ and the applied voltage for the P to AP and P to AP switching.
neuron. The numbers of on-centered receptive field array are weighted according to the Manhattan Distance [15] from the center of the window. Then, a 2D convolution operation was conducted and output a color image with analog values. The second step is to convert the color images into spike trains. Image data calculated from convolution operation is an analog value and can't be fed into the input neuron layer directly. Since the neuron's membrane potential could only be modulated by the spike stimulus, we employ a spike-train encoder built by interpolation function as an interface between image data from the physical world and SNNs.
The type of encoding adopted here is rate coding, which suggests that the information is carried by the firing rate of the neuron. Hence, spike trains generated with an average firing rate lies between 1-200 Hz corresponding to the input images have been shown in Figure S10.
The third step is the construction of spiking neurons. In this work, we choose a widely accepted simplified Leaky-Integrate-Fire (LIF) model as neuron research object. As illustrated in Figure S11a, a two-layer network composed of a set of input and output neurons were connected through synapses. Along with the arrival of encoded image information, the input neurons generate discrete events or spikes and send them to the output neurons. The spikes from Figure S10. Illustration of the transition process from the input images to spike trains. Firstly, the original images representing three hand-written numbers went through a convolution operation by the receptive field array and be converted into the images more suitable for human eye perception. Each box of the 20×20 grid is corresponding to one pixel in the input image. The receptive field array will act as a convolution kernel sliding from beginning to the end of the grid. Then, the post-treatment images were fed to an encoder, which could convert them into spike stimulus for input neurons as shown in the right three scatter diagrams. The x-axe is the neuron number while y-axe is the spike time. Every point in the graph represent a spike.
all the pre-neurons are altered as per the associated weights W1~n, and summed up as shown in Figure S11a. The output after summation alters the membrane potential ( mem ) of the postneuron in a typical Leaky-Integrate-Fire fashion. [16] When the potential crosses a threshold value (Vth), the post-neuron emits a spike and then enters into a refractory period during which no new input is received and the mem remains constant. Meanwhile, the first firing output neuron performs lateral inhibition on the rest of neurons. Figure S11b gives an example of simulation results to verify the correctness of our model design.
The fourth step is mapping the spin-torque memristor device into synapses. The key is to determine the resistance change δ of synapses after each epoch of stimulus, depending on the spike signal received by each synapse. This change should be consistent with the characteristic of the real devices, including its stochastic nature. The model introduced in Section S5 is used.
In this model, parameters of device such as , , + , + , − , − is extracted from the experimentally measured results on real devices. Then δ depends only on three parameters, i.e., the present resistance R, the amplitude V and the duration P of spike applied on the synapse, multiplied by a noise factor. As given in the main text ( Figure 5 c-e), to realize the STDP function, two sequences of ramped voltage pulses are used as pre-spike and a couple of short pulses with opposite polarity are used as post-spike. Since the amplitued of these two spikes are both lower than the critical current (see fig. 2b in the paper), neither of them can a b Figure S11. a) A representative model for the Leaky-Integrate-Fire (LIF) neuron used in the SNN simulation. The input spikes generated by input neurons are modulated by the synaptic weights W1~n and summed up together. The summation output alters the membrane potential of the post-neuron. The neuron will emit a spike if the membrane potential crosses a certain threshold. b) Simulation results for the 3 output neurons. Top panel shows the membrane potential of the winner neuron, exhibiting the typical LIF dynamics. No more spikes are generated during the refractory period until the device is reset to its initial potential. Middle and bottom panel shows the membrane potential of the other two output neurons that were inhibited by the fire of winner neuron.
change the resitance independently. Effectively, δ depends only on the peak of the across voltage V of the two spikes, as shown in the bottom panel in Figure 5c, which is determined by the delay time Δt between the pre-and post-spike. According to the waveform used in our experiments shown in Figure 5c, the relationship between V and Δt can be expressed as, = { 0.00187 × ∆ + 0.5525 ∆ < 0 0.00199 × ∆ − 0.5544 ∆ > 0 (S -12) where the unit of V is V and ∆ is μs. This relationship is imported in the SNN simulation. After combining the Equation S10, Equation S11 and Equation S12, the STDP of the memristor can be adequately described and integrated into the simulation.
With the behavior model of the spin-torque memristor, the STDP weight update rule and the SNN framework described above, we have conducted a simulation with unsupervised learning and the result is shown in Figure 6c.

Section 7. stability of the intermediate state at room temperature
The stability of intermediate state against the magnetic field at room temperature is measured and shown in figure S12. It can be seen that the field needed to destroy the intermediate state drops a little compared with that at low temperature ( fig.2d) because of thermal fluctuation. [17] SI References