Synchronisation of spin-orbit torque oscillators in disordered nanomagnet arrays

Abstract

Mutually synchronized spin-orbit-torque nano-oscillators are attractive microwave sources and candidate physical-computing elements, but array-scale coherence is limited by fabrication-induced frequency dispersion, nonlinear amplitude-phase coupling, thermal phase noise, and nonuniform coupling. We report a theory and simulation study of 10 x 10 arrays of elliptical CoFeB nanomagnets driven by damping-like spin Hall torque from a Pt underlayer. Stochastic Landau-Lifshitz-Gilbert simulations were combined with a calibrated nonlinear auto-oscillator phase model to quantify how anisotropy, dimension, damping, and spin Hall angle disorder affect locking thresholds, linewidth narrowing, cluster formation, and phase-slip statistics. In the clean array, nearest-neighbour spin-wave and dipolar coupling produced a global Kuramoto order parameter R = 0.94 at a current density of 7.4 x 10^11 A m^-2 and reduced the simulated microwave linewidth from 71 MHz for an isolated oscillator to 8.6 MHz for the array output. With 4% anisotropy disorder and 3% lateral-size disorder, global locking was replaced by two to four coherent clusters unless a weak common electrical feedback path was added. The reduced phase model reproduced the micromagnetic locking boundary to within 6% in current density, but under-predicted rare thermally activated phase slips at cluster edges. The results indicate that disorder tolerance in spin-orbit oscillator arrays depends less on the mean coupling strength than on the overlap between the coupling graph, the nonlinear frequency shift, and the spatial correlation length of fabrication variation.

Introduction

Spin-transfer and spin-orbit torques provide a route to nanoscale microwave oscillators whose frequency, phase, and amplitude can be controlled electrically. The basic spin-transfer mechanism was established in the theoretical work of Slonczewski and Berger [1,2], and microwave oscillations in nanomagnets and point contacts were later observed experimentally in spin-polarized current devices [3,4]. A central attraction of these devices is that their small volume and nonlinear frequency tunability make them easy to bias, tune, and couple. The same properties also make them sensitive to disorder, temperature, and load impedance.

Synchronization has been part of the spin-torque oscillator problem from the beginning. Mutual phase locking of two microwave spin-torque nano-oscillators was demonstrated in closely spaced point-contact devices [5,6], and nonlinear auto-oscillator theory subsequently provided a compact language for describing phase locking, linewidth, amplitude relaxation, and frequency pulling [7,8,9]. For large arrays, the analogy with phase-oscillator networks and the Kuramoto model is useful but incomplete: spintronic oscillators have strong amplitude-phase coupling, current-dependent nonlinearity, magnetostatic fields, spin-wave propagation, and electrical feedback paths that can all shift the effective coupling phase [10].

Spin-orbit torque broadens this design space. Heavy-metal/ferromagnet bilayers can generate damping-like torque through the spin Hall effect [11], enabling nano-oscillators driven by pure spin current rather than a charge current through the magnetic free layer [12]. Spin Hall nano-oscillators have been synchronized to external microwave signals [13], mutually synchronized over long distances [14], and organized into two-dimensional arrays with neuromorphic-computing motivation [15]. At the same time, modelling work has shown that oscillator arrays can host cluster states, vortices, and topological phase textures rather than one simple in-phase mode [16].

The present article is a fictional but physically bounded simulation study of a specific question: how much fabrication disorder can a two-dimensional spin-orbit-torque oscillator array tolerate before global synchronization fails? We focus on elliptical nanomagnets coupled by dipolar fields, weak spin-wave leakage through a shared ferromagnetic film edge, and optional electrical feedback through a common microwave load. The goal is not to claim an experimental device yield. It is to identify which disorder channels matter most and where reduced phase models remain trustworthy.

Array geometry and magnetic model

The simulated device consists of 100 elliptical CoFeB nanomagnets, arranged as a 10 x 10 square array on a 5 nm Pt underlayer. Each nominal magnet is 120 nm x 60 nm x 1.6 nm, with centre-to-centre pitch of 260 nm. The long axes are aligned along x, and the equilibrium magnetisation lies close to the in-plane easy axis at zero current. The magnetic parameters are saturation magnetisation Ms = 1.05 MA m^-1, exchange stiffness A = 15 pJ m^-1, Gilbert damping alpha = 0.012, interfacial perpendicular anisotropy chosen to give an effective in-plane anisotropy field of 38 mT, and spin Hall angle theta_SH = 0.085. These values are not fitted to a single published stack; they are intended to represent a plausible Pt/CoFeB-type oscillator design.

Magnetisation dynamics were computed with a stochastic Landau-Lifshitz-Gilbert equation including damping-like spin Hall torque, a smaller field-like torque equal to 8% of the damping-like term, Oersted field from the Pt strip, dipolar coupling between all magnets, and thermal noise at 300 K. The computational cell size was 4 x 4 x 1.6 nm^3, and the integration time step was adaptive with a maximum step of 0.5 ps. Each simulation ran for 450 ns after a 50 ns turn-on transient. The microwave signal was estimated from the array-averaged anisotropic magnetoresistance proxy m_x m_y and Fourier transformed after applying a Hann window.

Fabrication disorder was represented by independent and spatially correlated perturbations to four parameters: lateral dimensions, effective anisotropy, Gilbert damping, and spin Hall angle. Lateral dimensions were sampled by perturbing the major and minor axes before meshing. Anisotropy and spin Hall angle disorder were sampled from normal distributions with optional Gaussian spatial correlation length lambda_c. Unless otherwise stated, the baseline disordered case used standard deviations of 3% in lateral size, 4% in anisotropy, 2% in damping, and 3% in spin Hall angle, with lambda_c = 2 pitches. This produces free-running frequency spreads of 210-290 MHz near the operating current, comparable to the linewidth and detuning scales that control synchronization in earlier oscillator studies [5,7,8].

Reduced phase model

A reduced nonlinear auto-oscillator model was built from single-oscillator micromagnetic calibration. For each disorder realization, isolated magnets were simulated over current densities from 5.8 x 10^11 to 8.4 x 10^11 A m^-2. The resulting amplitude, frequency, amplitude-relaxation rate, and nonlinear frequency shift were fitted to the Slavin-Tiberkevich auto-oscillator form [8,9]. The fitted nonlinear frequency shift was negative over the full current range, with N/2pi between -0.8 and -1.4 GHz depending on anisotropy and aspect ratio. This nonlinearity is important because it converts amplitude perturbations caused by neighbours into phase shifts.

The phase dynamics were then written as dphi_i/dt = omega_i + sum_j K_ij sin(phi_j - phi_i - beta_ij) + eta_i(t), where omega_i is the calibrated natural frequency, K_ij is the coupling magnitude, beta_ij is a coupling phase lag, and eta_i is a thermal phase-noise term inferred from the isolated linewidth. Nearest-neighbour K_ij values were obtained from two-oscillator micromagnetic simulations by measuring injection-locking cones as a function of detuning. Next-nearest-neighbour and longer-range terms were scaled by the dipolar field plus an exponentially decaying spin-wave contribution. This is closer to a frustrated Kuramoto network than to an all-to-all phase model [10,16].

A second model variant included common electrical feedback. The array was connected to a simplified microwave load whose current modulation is proportional to the summed magnetoresistive signal with delay tau_e = 18 ps. The feedback strength was set below the value required to lock isolated oscillators by itself. Its role is to provide weak long-range coupling that supplements local magnetic interactions, similar in spirit to electrical coupling schemes proposed and demonstrated for spin-torque oscillator networks [17]. The electrical path is intentionally idealised: the model includes phase delay and load impedance but not a full distributed interconnect.

Simulation protocol and observables

For each disorder amplitude, 120 independent arrays were simulated with the reduced phase model. A subset of 24 arrays per disorder condition was then re-simulated micromagnetically to test the reduced model. Current density was swept upward in 0.1 x 10^11 A m^-2 increments, and the final 250 ns of each trajectory was used for analysis. Global coherence was quantified by the Kuramoto order parameter R = |N^-1 sum_j exp(i phi_j)|. We labelled a run globally locked when the time-averaged R exceeded 0.85 and the standard deviation of instantaneous oscillator frequencies was below 25 MHz. These thresholds are arbitrary but conservative relative to the clean-array linewidth.

We also tracked cluster order rather than treating synchronization as binary. Oscillators were assigned to the same cluster if their mean frequencies differed by less than 20 MHz and their phase difference remained bounded over the analysis window. Phase-slip events were counted when a nearest-neighbour phase difference crossed pi and did not immediately return within two oscillation periods. The output linewidth was extracted from the full-width at half-maximum of the principal spectral peak after subtracting the finite-window resolution determined from a deterministic clean-array run.

The principal sensitivity metric was the locking current J_lock, defined as the lowest current density at which at least 80% of disorder realizations satisfied the global-locking criterion. We report J_lock together with locked fraction, mean number of clusters, linewidth, and phase-slip rate. The distinction matters because a partially locked array can have useful microwave power and narrow local spectra even when its global order parameter is low.

Results

The clean array enters a globally synchronized state at J = 7.1 x 10^11 A m^-2 in the reduced model and at 7.2 x 10^11 A m^-2 in the full micromagnetic simulation. At J = 7.4 x 10^11 A m^-2, the micromagnetic array has R = 0.94 +/- 0.02, a mean oscillation frequency of 7.82 GHz, and a simulated output linewidth of 8.6 MHz. An isolated oscillator at the same current has a linewidth of 71 MHz under the same thermal-noise model. The linewidth reduction is smaller than ideal 1/N averaging because phase noise remains spatially correlated over several lattice spacings and because edge oscillators have a slightly different effective field.

Disorder changes the transition from sharp to gradual. With 2% anisotropy disorder and no size disorder, the locked fraction at 7.4 x 10^11 A m^-2 remains 0.91. With the baseline 4% anisotropy and 3% size disorder, it falls to 0.42 under local magnetic coupling alone. The typical state is not incoherent noise. It consists of two to four frequency-locked clusters separated by rows or corners where natural-frequency detuning exceeds the local coupling bandwidth. Cluster boundaries are preferentially located where correlated anisotropy gradients cross the array diagonal, demonstrating that spatial correlation in fabrication variation matters even when the one-point disorder distribution is unchanged.

Among the four disorder channels, effective anisotropy is dominant. Holding the root-mean-square free-running frequency spread fixed, anisotropy disorder produces more phase slips than spin Hall angle disorder because it changes both natural frequency and nonlinear frequency shift. Size disorder is second most important because aspect-ratio changes alter demagnetising fields and edge-mode character. Damping disorder affects the locking threshold weakly unless it is correlated with spin Hall angle disorder, in which case some magnets sit close to the onset current and behave as noise sources. This hierarchy is consistent with nonlinear auto-oscillator theory: detuning alone does not determine locking when the amplitude relaxation and nonlinear frequency shift also vary [7,8,9].

The reduced phase model captures the main locking boundary. Across all disorder amplitudes tested, its predicted J_lock differs from the micromagnetic result by a mean of 5.8%. The agreement is best for uncorrelated disorder and for arrays in which all oscillators remain on the same orbit type. It is worse near cluster boundaries, where edge oscillators intermittently excite nonuniform modes that are not represented by a single phase variable. In those cases the phase model under-predicts the phase-slip rate by a factor of 1.6-2.3, even though it still predicts the mean cluster count correctly.

Adding weak common electrical feedback substantially improves disorder tolerance. In the baseline disorder case, the globally locked fraction at 7.4 x 10^11 A m^-2 increases from 0.42 to 0.76, and the median linewidth narrows from 27 MHz to 13 MHz. The feedback does not simply force all oscillators into phase. Depending on the delay, it selects either an almost in-phase mode or a shallow phase-gradient mode with opposite edges shifted by 0.3-0.5 rad. Delays above 45 ps destabilise the in-phase branch and produce breathing clusters. This sensitivity to feedback phase agrees qualitatively with the broader oscillator-network literature and with long-range coupling observations in spin Hall and spin-torque arrays [14,15,17].

Discussion

The main result is that array coherence is controlled by the competition between disorder-induced detuning and the structure of the coupling network, not by a single average coupling value. A design with strong local coupling can still fail if correlated fabrication variation creates a smooth frequency gradient across the lattice. Conversely, weak long-range electrical feedback can rescue global order when it bridges cluster boundaries. This explains why the same disorder amplitude can produce either a globally locked array or a multi-cluster state depending on its spatial correlation length.

The results also clarify when a Kuramoto-like model is useful. Once the single-oscillator orbit and coupling phase are calibrated, the phase model is accurate enough to screen hundreds of disorder realizations and identify likely locking thresholds. It fails most visibly when oscillators switch between spatial modes or when thermal fluctuations drive rare slips at cluster edges. This limitation is not surprising: the classical Kuramoto framework was developed for phase oscillators [10], while spintronic auto-oscillators remain amplitude-active, nonlinear magnetic objects. The reduced model should therefore be treated as a design accelerator, not as a replacement for selected micromagnetic validation.

For device design, the practical message is that disorder budgets should be expressed in frequency and nonlinear-parameter space, not only in lithographic dimensions. A 3% width variation is tolerable in some stacks and severe in others depending on the slope df/dw and on how the nonlinear frequency shift changes with aspect ratio. Similarly, spin Hall angle nonuniformity matters most when it pushes a subset of oscillators close to threshold. Reporting only mean current density and mean frequency would hide these failure modes.

Several limitations are important. The simulated magnets are simpler than real nanoconstriction spin Hall oscillators, whose mode profiles can be strongly nonuniform and current-confined. Joule heating is represented only through the 300 K thermal-noise term and not through spatial temperature gradients or temperature-dependent material parameters. The common electrical feedback model is a lumped approximation and does not include parasitic capacitance, packaging modes, or impedance mismatch. The study also assumes that the disorder distributions are known. In an experimental array, they would need to be inferred from microscopy, single-device calibration, or Bayesian fitting to spectra. These limitations mean that the numerical thresholds should not be read as process specifications, but the qualitative hierarchy of disorder channels is likely robust.

Conclusion

We simulated synchronization in two-dimensional spin-orbit-torque oscillator arrays under realistic but fictional fabrication disorder. Clean 10 x 10 arrays locked globally and showed substantial linewidth narrowing, but moderate anisotropy and size disorder produced coherent clusters rather than uniform phase locking. A calibrated nonlinear phase model reproduced the micromagnetic locking boundary within 6% in current density, while missing some rare phase slips associated with mode changes at cluster boundaries.

The analysis suggests three design rules. First, anisotropy and aspect-ratio control deserve higher priority than damping uniformity once all magnets are comfortably above onset. Second, spatial correlation in disorder should be measured because gradients are more damaging than uncorrelated noise with the same variance. Third, weak long-range electrical feedback can improve disorder tolerance, but only if its phase delay is included in the design. These conclusions align with the broader view of spin-torque and spin Hall oscillators as nonlinear building blocks whose array behaviour emerges from both magnetic dynamics and circuit-level coupling [18].

Data and code availability

The supplementary archive contains all micromagnetic parameter files, reduced-model scripts, disorder seeds, coupling matrices, and spectral-analysis notebooks used for this study. The reduced phase-model code was tested with Python 3.10, NumPy 1.24, SciPy 1.10, and h5py 3.8. Full micromagnetic trajectories are provided for the clean array and for representative locked, clustered, and unlocked disordered arrays; the complete ensemble is represented by compressed phase and amplitude time series because of file size.