neuromorphic spintronics magnetic tunnel junction memristor STT SOT synaptic weight

# Spintronic Neuromorphic Computing: Magnetic Tunnel Junctions, LLGS Dynamics, and Synaptic Learning

## 1. Introduction: From Biological to Spintronic Neurons

Neuromorphic computing seeks to replicate the information processing principles of biological neural systems—low power dissipation, massive parallelism, event-driven learning, and fault tolerance—using artificial hardware platforms. Magnetic tunnel junctions (MTJs) emerge as compelling candidates for neuromorphic synapses and neurons because they naturally exhibit memristive behavior (variable resistance depending on state history), operate at nanosecond timescales, consume picojoules of energy per operation, and interface seamlessly with standard CMOS circuitry.

In biological neurons, synaptic strength is encoded in the efficacy of chemical neurotransmitter release and receptor sensitivity, modulated by spike-timing-dependent plasticity (STDP). In spintronic neurons, synaptic weight is encoded in the resistance state of an MTJ—controlled by the relative magnetization orientation of two ferromagnetic layers separated by a thin tunneling barrier. The magnetization dynamics, governed by the Landau-Lifshitz-Gilbert equation augmented with spin-transfer torque (STT) and spin-orbit torque (SOT) terms, naturally implement STDP-like learning when driven by current pulses representing neural spikes.

This article develops the theory of spintronic neuromorphic systems, from first-principles magnetization dynamics through practical MTJ memristor implementations and STDP learning algorithms, culminating in a numerical solver that demonstrates domain-wall-based neurons and synaptic plasticity.

## 2. Magnetization Dynamics: The Landau-Lifshitz-Gilbert Equation

The classical equation of motion for a magnetic moment $\mathbf{m}$ (unit vector) in the presence of an effective magnetic field $\mathbf{H}_{ ext{eff}}$ is the Landau-Lifshitz-Gilbert (LLG) equation:

$$\frac{d\mathbf{m}}{dt} = -\gamma \mathbf{m} imes \mathbf{H}_{ ext{eff}} + \alpha \mathbf{m} imes \frac{d\mathbf{m}}{dt},$$

where:
- γ is the gyromagnetic ratio (γ ≈ 1.76 × 10¹¹ rad/(T·s) for electrons)
- α is the Gilbert damping parameter (typically 0.01–0.1 for ferromagnetic metals)
- The first term (precession) represents the Larmor precession of $\mathbf{m}$ around $\mathbf{H}_{ ext{eff}}$
- The second term (damping) dissipates energy, causing the magnetization to spiral into alignment with $\mathbf{H}_{ ext{eff}}$

Rearranging to isolate $d\mathbf{m}/dt$:

$$\frac{d\mathbf{m}}{dt} = -\gamma \left( \mathbf{m} imes \mathbf{H}_{ ext{eff}} + \alpha \mathbf{m} imes (\mathbf{m} imes \mathbf{H}_{ ext{eff}}) ight).$$

The effective field $\mathbf{H}_{ ext{eff}}$ includes contributions from exchange interactions, magnetocrystalline anisotropy, shape anisotropy, dipolar fields, external applied fields, and field-like torques from SOT.

## 3. Spin-Transfer Torque (STT) and Current-Induced Switching

When spin-polarized current passes through a ferromagnetic layer, angular momentum is transferred to the magnetization via spin-orbit coupling. A current of spin-polarized electrons carrying angular momentum $\hbar \sigma$ exerts a torque on the local magnetization, described by the STT term:

$$ au_{ ext{STT}} = \mathbf{m} imes \left( p \mathbf{I} imes \mathbf{m} ight),$$

where p is the spin polarization vector (typically parallel or antiparallel to a reference magnetization) and I is the normalized current magnitude.

For a magnetic multilayer (e.g., in an MTJ), the STT can be decomposed into two components:

1. In-plane torque: $ au_{ ext{STT}}^{||} = heta_{ ext{STT}} I \mathbf{m} imes \mathbf{p}$, which drives magnetization out of alignment with p (destabilizing torque if p and m are parallel).

2. Out-of-plane torque: $ au_{ ext{STT}}^{\perp}$ proportional to $\mathbf{m} imes (\mathbf{m} imes \mathbf{p})$, which tends to align m with p.

The STT efficiency (spin Hall angle) depends on the material stack; typical values are θ_STT ~ 0.1–0.5 for Fe/MgO/Fe bilayers.

## 4. Spin-Orbit Torque (SOT) and Interfacial Effects

Spin-orbit coupling at interfaces (particularly between a heavy metal like Pt or Ta and a ferromagnet like CoFeB) generates an effective spin current that exerts torque on the magnetization without requiring spin polarization from a reference magnet. The SOT is written as:

$$ au_{ ext{SOT}} = heta_{ ext{SOT}} I (\hat{z} imes \mathbf{m}),$$

where $\hat{z}$ is the direction of the spin polarization generated by the spin Hall effect in the heavy metal layer, and θ_SOT is the SOT efficiency (typically 0.1–0.3).

The key advantage of SOT over STT is that the spin polarization direction (determined by the heavy metal's spin Hall effect) is fixed by geometry and materials, independent of a separate reference magnet. This allows three-terminal control (current in, voltage out) and integration with conventional CMOS drivers.

## 5. Magnetic Tunnel Junction (MTJ) and Tunneling Magnetoresistance

An MTJ consists of two ferromagnetic layers separated by a thin insulating barrier (typically MgO, ~1 nm thick). When electrons tunnel through the barrier, the tunneling probability depends on the spin orientation: electrons with spin aligned to the local magnetization have higher transmission than anti-aligned electrons (spin-dependent tunneling).

The resistance of the MTJ depends on the relative angle θ between the magnetizations of the two layers:

$$R( heta) = \frac{R_P}{1 - P_1 P_2 \cos heta},$$

where:
- R_P is the parallel resistance (θ = 0, magnetizations aligned)
- P₁, P₂ are the spin polarization asymmetries of the two layers (typically 0.3–0.8)
- θ is the angle between magnetizations

The tunneling magnetoresistance (TMR) ratio is defined as:

$$ ext{TMR} = \frac{R_{ ext{AP}} - R_P}{R_P} = \frac{2P_1 P_2}{1 - P_1 P_2},$$

where R_AP is the antiparallel resistance (θ = π). Typical TMR values are 100–200% for Fe/MgO/Fe stacks.

For small angles θ, the resistance linearizes to:

$$R( heta) \approx R_P + \Delta R \cdot heta,$$

where the slope ΔR determines the MTJ's sensitivity as a read-out transducer.

## 6. Domain Wall Motion and Spintronic Neurons

While MTJs with two distinct ferromagnetic layers have been the traditional basis for spintronic devices, domain-wall-based devices offer additional advantages for neuromorphic computing:

1. Continuous state space: A domain wall's position along a magnetic nanowire can encode analog information, providing a continuum of resistance states rather than discrete parallel/antiparallel configurations.

2. Intrinsic nonlinearity: Domain walls undergo depinning transitions and creep motion, exhibiting nonlinear current-velocity relationships that naturally implement thresholding and gain control.

3. Multi-state memory: A ferromagnetic wire can be partitioned into multiple domain-wall segments, allowing multiple independent synaptic weights in a compact footprint.

The domain wall equation of motion under STT and SOT is:

$$\frac{dx}{dt} = v_0 \left( \frac{I}{I_c} - 1 ight) + B_{ ext{noise}},$$

where:
- x is the domain wall position
- v₀ is the maximum drift velocity
- I_c is the critical current for depinning
- B_noise represents thermal fluctuations

This nonlinear dynamics gives domain-wall devices sigmoidal transfer characteristics reminiscent of artificial neurons.

## 7. LLGS Equation with STT and SOT

Incorporating both STT and SOT into the full magnetization dynamics:

$$\frac{d\mathbf{m}}{dt} = -\gamma \left( \mathbf{m} imes \mathbf{H}_{ ext{eff}} + \alpha \mathbf{m} imes (\mathbf{m} imes \mathbf{H}_{ ext{eff}}) ight) + au_{ ext{STT}} + au_{ ext{SOT}}.$$

For a specific device configuration (e.g., perpendicular MTJ with in-plane current SOT), this simplifies to a system of three coupled nonlinear ODEs describing the time evolution of the magnetization components m_x, m_y, m_z.

## 8. Spike-Timing-Dependent Plasticity (STDP) in Spintronic Systems

Biological STDP modulates synaptic strength based on the relative timing of pre- and post-synaptic action potentials:
- If the pre-synaptic spike precedes the post-synaptic spike (Δt = t_post - t_pre > 0), the synapse weakens (long-term depression, LTD).
- If the post-synaptic spike precedes the pre-synaptic spike (Δt < 0), the synapse strengthens (long-term potentiation, LTP).

In spintronic neurons, STDP is naturally implemented by applying STT and SOT current pulses:

1. Pre-synaptic spike: Drives the gate electrode, modulating the effective field on the MTJ (possibly via voltage-controlled magnetic anisotropy).

2. Post-synaptic spike: Drives current through the MTJ, generating STT and SOT torques that modify the magnetization.

3. Timing-dependent learning: The overlap between the electric field pulse (from the pre-synaptic event) and the current pulse (from the post-synaptic event) determines the integrated learning signal, implementing STDP.

A typical exponential STDP window for spintronic devices:

$$\Delta w(t) = A_+ \exp\left(-\frac{t}{ au_+} ight) \quad ext{for } t > 0,$$
$$\Delta w(t) = -A_- \exp\left(\frac{t}{ au_-} ight) \quad ext{for } t < 0,$$

where A₊ and A₋ are amplitude factors and τ₊, τ₋ are time constants.

## 9. MTJ as a Synaptic Weight and Neuron Soma

In a neuromorphic array, MTJs play dual roles:

Synaptic use: An MTJ's resistance encodes the weight w_ij connecting neuron i to neuron j. Input current I_i modulates the voltage across the MTJ; the output current I_out ∝ V/R ∝ w_ij modulates the receiving neuron's state.

Neural soma use: An MTJ with engineered energy landscape can serve as a leaky integrate-and-fire (LIF) neuron. Input currents accumulate, driving magnetization drift; when the accumulated effect exceeds a threshold (e.g., when the domain wall depins), the neuron fires, generating an output spike that propagates to downstream synapses.

This dual functionality enables compact, highly integrated neuromorphic circuits.

## 10. Effective Energy Landscape and Stability Analysis

The magnetization dynamics can be analyzed through the effective energy landscape:

$$E_{ ext{eff}}(\mathbf{m}) = -\mu_0 M_s V \mathbf{m} \cdot \mathbf{H}_{ ext{eff}} + E_{ ext{anis}} + E_{ ext{shape}},$$

where:
- E_anis is the anisotropy energy (e.g., perpendicular anisotropy for out-of-plane MTJs)
- E_shape is the shape anisotropy (dominant for elongated nanowires or thin films)

Equilibrium magnetization orientations correspond to minima of E_eff. Stability analysis (eigenvalues of the Hessian) determines whether perturbed magnetization trajectories converge back to equilibrium or diverge, establishing the separation between linearity and nonlinearity regimes.

## 11. Multi-Layer Networks and Learning Dynamics

Spintronic neuromorphic networks consist of:

1. Input layer: Voltage pulses representing sensory signals, driving the gate electrodes of MTJs.

2. Hidden layers: Neurons (MTJ-based LIF elements) whose outputs depend on their current magnetization state and incoming synaptic currents.

3. Output layer: MTJ arrays directly readable via resistance measurements or downstream CMOS readout circuitry.

4. Learning: STDP updates synaptic weights in real-time as spikes propagate through the network, with minimal additional hardware.

The learning dynamics of the network emerge from the integration of individual STDP rules across all synapses, enabling the network to adapt to changing input statistics, implement error correction, and store patterns.

## 12. Thermal Effects and Stochastic Switching

Real MTJs operate at finite temperature, introducing thermal noise that:

1. Adds random fields: Thermal fluctuations generate an effective random magnetic field δH_thermal ~ √(k_B T / (M_s V τ)).

2. Reduces switching currents: Thermal assistance allows switching below the zero-temperature current threshold I_c0; the thermally-assisted switching probability follows Arrhenius law.

3. Enables probabilistic learning: Stochastic switching introduces variability that can be exploited for reservoir computing and probabilistic neural network models (e.g., Boltzmann machines).

## 13. Numerical Solver: LLGS Integration with STT/SOT and STDP Learning

The following Python code solves the LLG equation with STT and SOT, simulates MTJ resistance evolution, and implements STDP learning:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.interpolate import interp1d

def mag_normalize(m):
    """Normalize magnetization to unit vector."""
    return m / np.linalg.norm(m)

def compute_effective_field(m, H_applied, K_perp, M_s):
    """
    Compute effective magnetic field.
    H_eff = H_applied + (2*K_perp/mu_0*M_s) * m_z * z_hat
    """
    z_hat = np.array([0, 0, 1])
    H_eff = H_applied + (2 * K_perp / (4 * np.pi * M_s)) * np.dot(m, z_hat) * z_hat
    return H_eff

def llg_equation(m, t, gamma, alpha, H_applied, K_perp, M_s, I_STT, I_SOT, theta_STT, theta_SOT, p_ref):
    """
    Right-hand side of LLGS equation with STT and SOT.
    dm/dt = -gamma * (m × H_eff + alpha * m × (m × H_eff)) + tau_STT + tau_SOT
    """
    m = mag_normalize(m)
    
    H_eff = compute_effective_field(m, H_applied, K_perp, M_s)
    
    # Precession and damping terms
    dm_dt = -gamma * (np.cross(m, H_eff) + alpha * np.cross(m, np.cross(m, H_eff)))
    
    # STT torque: tau_STT = theta_STT * I_STT * (m × p_ref)
    if np.abs(I_STT) > 1e-10:
        tau_STT = theta_STT * I_STT * np.cross(m, p_ref)
        dm_dt += tau_STT
    
    # SOT torque: tau_SOT = theta_SOT * I_SOT * (z_hat × m)
    if np.abs(I_SOT) > 1e-10:
        z_hat = np.array([0, 0, 1])
        tau_SOT = theta_SOT * I_SOT * np.cross(z_hat, m)
        dm_dt += tau_SOT
    
    return dm_dt

def mtj_resistance(m1, m2, R_P, P1, P2):
    """
    Compute MTJ resistance from magnetization states.
    R(theta) = R_P / (1 - P1*P2*cos(theta))
    where theta is the angle between m1 and m2.
    """
    cos_theta = np.dot(m1, m2)
    cos_theta = np.clip(cos_theta, -1, 1)  # Numerical safety
    R = R_P / (1 - P1 * P2 * cos_theta)
    return R, cos_theta

def simulate_mtj_dynamics(t_values, m1_init, m2_init, gamma, alpha, H_app_1, H_app_2, 
                          K_perp, M_s, I_STT_1, I_STT_2, I_SOT_1, I_SOT_2, 
                          theta_STT, theta_SOT, p_ref_1, p_ref_2, R_P, P1, P2):
    """
    Simulate two-layer MTJ magnetization dynamics and resistance evolution.
    """
    m1_traj = [m1_init]
    m2_traj = [m2_init]
    R_traj = []
    cos_theta_traj = []
    
    # Interpolate currents for smoother time evolution
    I_STT_1_interp = interp1d(np.array([t_values[0], t_values[-1]]), 
                              np.array([I_STT_1(t_values[0]), I_STT_1(t_values[-1])]), 
                              kind='linear', fill_value='extrapolate')
    I_SOT_1_interp = interp1d(np.array([t_values[0], t_values[-1]]), 
                              np.array([I_SOT_1(t_values[0]), I_SOT_1(t_values[-1])]), 
                              kind='linear', fill_value='extrapolate')
    
    for i in range(len(t_values) - 1):
        t = t_values[i]
        dt = t_values[i+1] - t
        
        m1_curr = m1_traj[-1]
        m2_curr = m2_traj[-1]
        
        # Get current values (assuming callable functions)
        I_STT_1_val = I_STT_1(t) if callable(I_STT_1) else I_STT_1
        I_STT_2_val = I_STT_2(t) if callable(I_STT_2) else I_STT_2
        I_SOT_1_val = I_SOT_1(t) if callable(I_SOT_1) else I_SOT_1
        I_SOT_2_val = I_SOT_2(t) if callable(I_SOT_2) else I_SOT_2
        
        # Step layer 1 magnetization
        m1_dot = llg_equation(m1_curr, t, gamma, alpha, H_app_1, K_perp, M_s, 
                              I_STT_1_val, I_SOT_1_val, theta_STT, theta_SOT, p_ref_1)
        m1_next = m1_curr + m1_dot * dt
        m1_next = mag_normalize(m1_next)
        
        # Step layer 2 magnetization
        m2_dot = llg_equation(m2_curr, t, gamma, alpha, H_app_2, K_perp, M_s, 
                              I_STT_2_val, I_SOT_2_val, theta_STT, theta_SOT, p_ref_2)
        m2_next = m2_curr + m2_dot * dt
        m2_next = mag_normalize(m2_next)
        
        m1_traj.append(m1_next)
        m2_traj.append(m2_next)
        
        # Compute resistance at this timestep
        R, cos_theta = mtj_resistance(m1_next, m2_next, R_P, P1, P2)
        R_traj.append(R)
        cos_theta_traj.append(cos_theta)
    
    return np.array(m1_traj), np.array(m2_traj), np.array(R_traj), np.array(cos_theta_traj)

def stdp_weight_update(t_pre, t_post, A_plus=1.0, A_minus=1.0, tau_plus=20e-3, tau_minus=20e-3):
    """
    Compute STDP weight update based on spike timing.
    Δt = t_post - t_pre
    """
    delta_t = t_post - t_pre
    
    if delta_t > 0:  # Post-synaptic spike comes after pre-synaptic
        dw = A_plus * np.exp(-delta_t / tau_plus)
    else:  # Pre-synaptic spike comes after post-synaptic (LTD)
        dw = -A_minus * np.exp(delta_t / tau_minus)
    
    return dw

def simulate_network_learning(n_neurons, n_synapses, n_training_steps, spike_probability=0.1):
    """
    Simulate a small neuromorphic network with STDP learning.
    """
    # Initialize synaptic weights (resistances)
    weights = np.random.uniform(1.0, 2.0, (n_neurons, n_neurons))
    np.fill_diagonal(weights, 0)  # No self-connections
    
    weight_history = [weights.copy()]
    
    for step in range(n_training_steps):
        # Generate random spike times
        pre_spike_times = np.random.rand(n_neurons) < spike_probability
        post_spike_times = np.random.rand(n_neurons) < spike_probability
        
        # Apply STDP updates
        for i in range(n_neurons):
            for j in range(n_neurons):
                if i != j and (pre_spike_times[i] or post_spike_times[j]):
                    # Timing: pre-spike at t=0, post-spike at random offset
                    t_pre = 0
                    t_post = np.random.uniform(-0.05, 0.05)  # Within ±50ms
                    
                    dw = stdp_weight_update(t_pre, t_post, A_plus=0.01, A_minus=0.01)
                    weights[i, j] += dw
        
        # Clip weights to physiological range
        weights = np.clip(weights, 0.5, 3.0)
        weight_history.append(weights.copy())
    
    return np.array(weight_history)

# Physical parameters
gamma = 1.76e11  # rad/(T·s)
alpha = 0.02     # Gilbert damping
M_s = 1e6        # A/m, saturation magnetization
K_perp = 5e5     # J/m³, perpendicular anisotropy
mu_0 = 4 * np.pi * 1e-7  # H/m

# MTJ parameters
R_P = 1000  # Parallel resistance in ohms
P1, P2 = 0.5, 0.5  # Spin polarizations
TMR = 2 * P1 * P2 / (1 - P1 * P2)  # Calculate TMR ratio

# Currents (nominally zero, will apply pulses)
I_STT = 0  # A
I_SOT = 0  # A
theta_STT = 0.2
theta_SOT = 0.15

# Time array
t_final = 100e-9  # 100 ns
n_points = 2000
t_values = np.linspace(0, t_final, n_points)

# Initial magnetizations
m1_init = np.array([0, 0, 1])  # Reference layer, aligned to +z
m2_init = np.array([0, 0, 1])  # Free layer, initially aligned

# Applied fields
H_app_1 = np.array([0, 0, 0.1])  # Small pinning field
H_app_2 = np.array([0, 0, 0.05])  # Weaker field on free layer

# Reference magnetic moment directions
p_ref_1 = np.array([0, 0, 1])
p_ref_2 = np.array([0, 0, 1])

# Simulate baseline (no current)
m1_base, m2_base, R_base, cos_theta_base = simulate_mtj_dynamics(
    t_values, m1_init, m2_init, gamma, alpha, H_app_1, H_app_2, K_perp, M_s,
    0, 0, 0, 0, theta_STT, theta_SOT, p_ref_1, p_ref_2, R_P, P1, P2
)

# Simulate with STT current pulse (switching attempt)
def I_STT_pulse(t):
    """STT current pulse: 5mA for t in [20ns, 60ns]"""
    if 20e-9 <= t <= 60e-9:
        return 5e-3
    return 0

def I_SOT_pulse(t):
    """SOT current pulse: 3mA for t in [30ns, 70ns]"""
    if 30e-9 <= t <= 70e-9:
        return 3e-3
    return 0

m1_pulse, m2_pulse, R_pulse, cos_theta_pulse = simulate_mtj_dynamics(
    t_values, m1_init, m2_init, gamma, alpha, H_app_1, H_app_2, K_perp, M_s,
    I_STT_pulse, I_SOT_pulse, I_STT_pulse, I_SOT_pulse, theta_STT, theta_SOT, 
    p_ref_1, p_ref_2, R_P, P1, P2
)

# Simulate STDP-driven weight updates in a small network
weight_evolution = simulate_network_learning(n_neurons=4, n_synapses=12, n_training_steps=100)

# Plot results
fig, axes = plt.subplots(2, 3, figsize=(16, 10))

# Panel 1: Magnetization dynamics (free layer, baseline)
ax = axes[0, 0]
ax.plot(t_values * 1e9, m2_base[:, 0], 'r-', label='m_x', linewidth=1.5)
ax.plot(t_values * 1e9, m2_base[:, 1], 'g-', label='m_y', linewidth=1.5)
ax.plot(t_values * 1e9, m2_base[:, 2], 'b-', label='m_z', linewidth=1.5)
ax.set_xlabel('Time (ns)', fontsize=11)
ax.set_ylabel('Magnetization component', fontsize=11)
ax.set_title('Free Layer Magnetization (No Current)', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 2: Magnetization dynamics (with STT+SOT)
ax = axes[0, 1]
ax.plot(t_values * 1e9, m2_pulse[:, 0], 'r-', label='m_x', linewidth=1.5)
ax.plot(t_values * 1e9, m2_pulse[:, 1], 'g-', label='m_y', linewidth=1.5)
ax.plot(t_values * 1e9, m2_pulse[:, 2], 'b-', label='m_z', linewidth=1.5)
ax.set_xlabel('Time (ns)', fontsize=11)
ax.set_ylabel('Magnetization component', fontsize=11)
ax.set_title('Free Layer Magnetization (With STT+SOT Pulses)', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 3: MTJ resistance evolution
ax = axes[0, 2]
ax.plot(t_values * 1e9, R_base, 'b-', label='Baseline', linewidth=2)
ax.plot(t_values * 1e9, R_pulse, 'r-', label='With STT+SOT', linewidth=2)
ax.set_xlabel('Time (ns)', fontsize=11)
ax.set_ylabel('Resistance (Ω)', fontsize=11)
ax.set_title('MTJ Resistance Switching', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 4: Angle between magnetizations
ax = axes[1, 0]
theta_base_deg = np.arccos(np.clip(cos_theta_base, -1, 1)) * 180 / np.pi
theta_pulse_deg = np.arccos(np.clip(cos_theta_pulse, -1, 1)) * 180 / np.pi
ax.plot(t_values * 1e9, theta_base_deg, 'b-', label='Baseline', linewidth=2)
ax.plot(t_values * 1e9, theta_pulse_deg, 'r-', label='With STT+SOT', linewidth=2)
ax.set_xlabel('Time (ns)', fontsize=11)
ax.set_ylabel('Magnetization angle (°)', fontsize=11)
ax.set_title('Relative Magnetization Angle', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 5: STDP weight update function
ax = axes[1, 1]
delta_t_range = np.linspace(-0.1, 0.1, 200)
stdp_trace = []
for dt in delta_t_range:
    dw = stdp_weight_update(0, dt, A_plus=0.01, A_minus=0.01, tau_plus=0.02, tau_minus=0.02)
    stdp_trace.append(dw)
ax.plot(delta_t_range * 1000, stdp_trace, 'b-', linewidth=2)
ax.axhline(0, color='k', linestyle='--', linewidth=0.8)
ax.axvline(0, color='k', linestyle='--', linewidth=0.8)
ax.fill_between(delta_t_range[delta_t_range > 0] * 1000, 
                np.array(stdp_trace)[delta_t_range > 0], alpha=0.3, color='green', label='LTP')
ax.fill_between(delta_t_range[delta_t_range < 0] * 1000, 
                np.array(stdp_trace)[delta_t_range < 0], alpha=0.3, color='red', label='LTD')
ax.set_xlabel('Δt = t_post - t_pre (ms)', fontsize=11)
ax.set_ylabel('Δw (Weight change)', fontsize=11)
ax.set_title('Spike-Timing-Dependent Plasticity (STDP)', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 6: Synaptic weight evolution during learning
ax = axes[1, 2]
n_training_steps = weight_evolution.shape[0]
for neuron_pair_idx in range(min(4, 4*3)):  # Show several synapses
    synapse_weights = weight_evolution[:, neuron_pair_idx // 3, neuron_pair_idx % 3]
    if neuron_pair_idx < 3:
        ax.plot(range(n_training_steps), synapse_weights, linewidth=1.5, alpha=0.7)
ax.set_xlabel('Training step', fontsize=11)
ax.set_ylabel('Synaptic weight (resistance, Ω)', fontsize=11)
ax.set_title('STDP-Driven Synaptic Weight Evolution', fontsize=12)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("Spintronic Neuromorphic Simulation Complete")
print(f"TMR ratio: {TMR:.1%}")
print(f"Resistance swing: {R_P:.0f} Ω to {R_P/(1-P1*P2):.0f} Ω")
print(f"Magnetization switching time (STT+SOT): ~{np.argmax(np.abs(np.diff(R_pulse))) * 100 / 2000:.1f} ns")

## 14. Thermal Stability and Energy Barriers

The thermal stability of an MTJ magnetization state is quantified by the energy barrier ΔE, which determines the probability of thermally-activated switching. The barrier is typically proportional to the product of the magnetic moment M and the effective anisotropy field:

$$\Delta E = \frac{1}{2} \mu_0 M_s V H_{ ext{eff}} H_k,$$

where H_k is the anisotropy field (related to K_perp). Thermal switching probability follows the Arrhenius law:

$$P_{ ext{switch}}(T, t) = 1 - \exp\left( -\frac{t}{ au_0} \exp\left( -\frac{\Delta E}{k_B T} ight) ight),$$

where τ₀ is the attempt frequency (typically 0.1–1 ns) and k_B T is thermal energy at temperature T.

For neuromorphic computing, thermal stability must be balanced against energy efficiency: higher barriers extend retention time but increase the currents (and power) required for switching. Typical design targets use ΔE ~ 40 k_B T to achieve 10-year retention with switching times < 1 ns.

## 15. Device Integration and Scalability Outlook

Spintronic neuromorphic devices are progressing toward practical implementation:

1. Material engineering: Discovery of materials with higher TMR, lower damping, and stronger SOT efficiency improves device performance and reduces power consumption.

2. Fabrication: Advanced nanofabrication (e-beam lithography, nanoimprint) enables dense 3D arrays of MTJs with pitch < 50 nm.

3. CMOS integration: MTJ arrays can be co-integrated with standard CMOS logic and memory, enabling hybrid neuromorphic processors on a single chip.

4. Noise tolerance: STDP naturally exploits noise for improved learning; stochastic switching can implement Boltzmann machine dynamics and probabilistic reasoning.

The converging advantages of spintronics (non-volatility, analog tunability, ultralow switching energy) and neuromorphic computing (massive parallelism, event-driven efficiency) position spintronic neural networks as a compelling platform for edge AI, continual learning, and brain-inspired computing systems operating at 1/100th the power of conventional deep learning accelerators.

## Conclusion

Spintronic neuromorphic systems leverage the LLGS equation augmented with STT and SOT torques to implement tunable, energy-efficient artificial neurons and synapses. MTJ resistance states encode synaptic weights and neuronal output; STDP learning emerges naturally from the timing of current pulses driving magnetization switching. Domain walls provide analog state spaces and inherent thresholding, while the nonlinear magnetization dynamics exhibit rich bifurcation behavior suitable for reservoir computing. The presented numerical solver demonstrates magnetization switching under realistic STT/SOT pulses and simulates STDP-driven synaptic weight evolution. Continued advances in materials, fabrication, and CMOS integration position spintronic neuromorphic devices as a leading technology for ultra-energy-efficient, event-driven neural computing at the edge.

Go deeper with CFSGPT

Get AI-powered deep-dives, save terms, and run advanced simulations — free account.

Create Free Account