Ultrafast Ising Machines using spin torque nano-oscillators Open Access

Combinatorial optimization problems are known for being particularly hard to solve on traditional von Neumann architectures. This has led to the development of Ising Machines (IMs) based on quantum annealers and optical and electronic oscillators, demonstrating speed-ups compared to central processing unit (CPU) and graphics processing unit (GPU) algorithms. Spin torque nano-oscillators (STNOs) have shown GHz operating frequency, nanoscale size, and nanosecond turn-on time, which would allow their use in ultrafast oscillator-based IMs. Here, we show using numerical simulations based on STNO auto-oscillator theory that STNOs exhibit fundamental characteristics needed to realize IMs, including in-phase/out-of-phase synchronization and second harmonic injection locking phase binarization. Furthermore, we demonstrate numerically that large STNO network IMs can solve Max-Cut problems on nanosecond timescales.

Combinatorial optimization (CO) problems appear in a wide range of scheduling, routing, and optimization situations in all of society. CO problems can be illustrated as the need for finding a global minimum in a multivariate energy landscape. As this is a very common situation in nature, solutions have been inspired by physical systems such as the Ising model.¹ The Hamiltonian associated with the Ising model can be written as

where N represents the number of spins, $Jj,j′$ is the coupling between spins j and $j′$⁠, s_j is the spin state +1 or −1 of spin j, and the last term describes a bias field h acting on all spins. Real-world CO problems can be mapped¹ onto the Hamiltonian in Eq. (1), which opens up the possibility of building a general purpose CO computing platform based on solving the Ising model.

Customized hardware architectures able to solve the Ising model Eq. (1) are called Ising Machines (IMs), and various implementations have been proposed. D-Wave quantum annealers based on Josephson junctions² have shown promising results but face similar challenges as general quantum computers: cryogenic temperatures requiring kWs of cooling power and extremely high production costs. IMs utilizing optoelectronic parametric oscillators address some challenges related to quantum systems, offering room temperature operation and all-to-all interconnections using optical feedback methods.^3–5 However, solutions to miniaturize these systems are challenging and are yet to be demonstrated. Alternatively, IMs inspired by simulated annealing and based on commercially available CMOS processes have shown promising results.^6–8 

A recent demonstration using LC oscillators and resistive crossbars proved that IM can be realized utilizing readily available electronic oscillators.^9,10 This approach is based on the similarities between the Ising Hamiltonian and the Lyapunov function describing the dynamics of a second harmonic injection locked (SHIL) network of coupled oscillators.

Networks of synchronized spintronic oscillators such as spin-torque nano-oscillators (STNOs) and spin-Hall nano-oscillators (SHNOs),¹¹ with characteristics including wide frequency tunability, nanoscale size, strongly non-linear behavior, and CMOS compatibility, have gained significant interest as building blocks for neuromorphic/unconventional computing.^12–19 These properties, along with the substantial effort focused on further enhancing the controllability of spintronic oscillator arrays,^20,21 make them interesting candidates for implementing IMs. It is, therefore, of great interest to explore whether networks of interacting STNOs, with their particularly strong and qualitatively different non-linearities, could be used as IMs. Here, we numerically explore this possibility.

The proposed IM consists of an array of coupled STNOs where SHIL is achieved by applying a small RF current, with twice the frequency of the free-running STNOs, in combination with the DC bias current to all the STNOs. To predict the behavior of such a network, we employ the universal auto-oscillator model proposed by Slavin and Tiberkevich,²² which describes the STNO dynamics in terms of the complex dimensionless spin wave amplitude c, where $p=|c2|$ is proportional to the STNO output power and $ϕ=arg(c)$ is the phase. The dynamical equations describing STNO j in an array under SHIL can be written as^22,23

or alternatively as two coupled real differential equations for the power p and phase $ϕ$ as

where $ωj(pj)=ω0+npj$ is the operating frequency, $ω0/2π≈4.2$ GHz is the ferromagnetic resonance frequency, $n/2π≈−3.44$ GHz is the nonlinear frequency shift coefficient, $Γ+,j(pj)=ΓG(1+Qpj)$ is the intrinsic magnetization damping (⁠$ΓG/2π≈252$ MHz), and $Q≈2.66$ is the nonlinear damping parameter. $Γ−,j(pj)=σIdc(1−pj)$ is the negative damping induced by the spin polarized direct current I_dc and $σ=εgμB/2eMsV$⁠, where $ε=0.7$ is the spin polarization, g is the spectroscopic Lande factor, μ_B is the Bohr magneton, e is the electron charge, $4πMs=9500$ Oe is the saturation magnetization, and $V=85×140×1.8$nm³ is the volume of the free layer. Parameters, including $n$ and ω₀, were calculated based on a CoFeB magnetic tunnel junction STNO,²⁴ as shown in Fig. 1(a), with an easy-plane magnetic anisotropy of H_A = 125 Oe and placed in an external in-plane magnetic field of H_ext = 300 Oe oriented at an angle of $ϕext=30°$ with respect to the anisotropy direction of the free layer. Moreover, the inter-layer coupling field is assumed as H_int = 125 Oe. All device parameters are taken from an analytical model developed in Ref. 25, which was successfully used to replicate the microwave signal properties of different independent magnetic tunnel junction STNO studies.^24,26

Injection locking is described using the frequency ω_e of the external RF current and its coupling strength, K_e to the STNOs. The last term on the right-hand side (RHS) describes the coupling of oscillator j to other oscillators $j′$ in the network through the coupling strength $Ωj,j′$ and phase $βj,j′$⁠. Common coupling mechanisms between STNOs include spin-wave,^{12–15,27,28} direct exchange,^16,19,29 dipolar,^30–32 and electrical^18,33–35 coupling, or even a combination of these mechanisms.^12,15,18,19 However, the type of coupling does not directly affect the model in Eq. (2) where only the strength and phase of the coupling mechanism are taken into account. Throughout the simulation results, we will focus on the behavior of the STNO system described by Eq. (2) without explicitly defining the type of coupling. However, a general purpose IM requires an all-to-all coupling, which limits the choice of physical coupling mechanisms.

Figure 1(b) presents simulation results for a single injection locked STNO [the last term on the RHS in Eq. (2) is neglected], where ω_e is swept around the STNO second harmonic, the so-called second harmonic injection locking (SHIL), with a fixed value of coupling to the external RF current, K_e = 50 MHz. For a large difference between $fe=ωe/2π$ and f_STNO, the STNO remains at the free running frequency of 3.2 GHz. However, when f_e approaches $2fSTNO$⁠, it starts perturbing the STNO, showing a typical frequency pulling between 2.9 GHz and 3.1 GHz, while the STNO completely locks to the external signal in the range of 3.1 GHz - 3.3 GHz. Similar to conventional electronic oscillators under SHIL, there are two stable locking states with a relative phase difference of $180°$⁠. This behavior is often referred to as phase binarization, which is critical for the realization of an oscillator-based IM. In the case of a single locked oscillator, the locking states $0°$/$180°$ depend purely on the initial conditions. However, for a network of coupled oscillators under SHIL, these two states can be used to represent spins in the Ising Hamiltonian as it has been previously demonstrated.⁹ 

Figures 1(c)–1(h) present the simulation results obtained by solving Eq. (2) for two identical coupled STNOs in the absence of SHIL [K_e = 0 in Eq. (2)]. STNO₁ is biased with a constant DC bias current resulting in $fSTNO1≈3.2$ GHz, while the DC bias/operating frequency of STNO₂ is swept. Figures 1(c) and 1(d) present the effect of increasing $Ωj,j′$⁠, which results in a wider frequency locking range. Moreover, the STNOs synchronize away from the mean of their individual operating frequencies as a consequence of the coupling phase $βj,j′$⁠.^22,35

The coupling phase characterizes the phase delay associated with the coupling mechanism, e.g., in the case of dipolar coupling $βj,j′=−π/2$⁠.²² Moreover, the coupling phase can be re-normalized to account for any other time delays of the coupling, such as, e.g., capacitive coupling.^36–38 A change in $βj,j′$ significantly modifies the locking bandwidth and the frequency at which the STNOs synchronize, as presented in Figure 1(d) and 1(e). Moreover, for the case of $βj,j′=−1.5π$ shown in Fig. 1(f), the STNOs synchronize in frequency but are out-of-phase (OOP) as presented in Fig. 1(h). This translates to a negative (antiferromagnetic) coupling $Ωj,j′$ in Eq. (2). For the conventional electronic oscillator, the OOP synchronization usually occurs in the fixed range of $0.5π<βj,j′<1.5π$⁠, while the nonlinearity in STNOs shifts this interval when Eq. (2) is solved.³⁹ Analytically, the effective coupling phase, $β̃j,j′$⁠, can be approximated as²² 

where ν is the normalized nonlinear frequency shift coefficient, calculated according to the macrospin theory.²² For the case presented in Fig. 1(h), $ν≈−2.64$⁠, which translates to an effective coupling phase of $β̃j,j′≈−1.1π$⁠, explaining why a robust OOP synchronization is achieved under these conditions. The combination of in-phase/OOP synchronization and SHIL phase binarization demonstrated above shows that STNOs exhibit all the desired characteristics to realize an IM.

Many NP-hard problems can be mapped onto the Ising model, but we will focus on a particular class of problems called Max-Cut, one of Karp's 21 NP-complete problems. Specifically, we focus on undirected, unweighted graphs with N vertices and E edges. The Max-Cut problem consists of partitioning the vertices into two subsets $s1/s2$ such that the number of edges crossing between s₁ and s₂ is as high as possible. This problem can be mapped to the Hamiltonian in Eq. (1) (assuming h = 0) where spins that settle to the states $↑/↓$ (+1/−1) belong to the subset $s1/s2$⁠. An optimal solution of a particular graph corresponds to the sets $s1/s2$ leading to the ground state of Eq. (1), which also represents the Max-Cut solution. An example graph, called a Möbius ladder graph with a size N = 16, is shown in Fig. 2(a) along with one possible ground state/Max-Cut solution.

In the oscillator-based IM, STNOs represent nodes/spins with the state +1/-1 determined by the phase $0°/180°$ and with graph edges programed into the couplings $Ωj,j′$⁠. Max-Cut problems require antiferromagnetic coupling $Ωj,j′$ between STNOs if nodes in a particular graph are connected and $Ωj,j′=0$ otherwise. The coupling phase is chosen as $βj,j′=−1.4π$ to yield an effective coupling phase of $β̃j,j′≈−π$⁠, which provides antiferromagnetic coupling in our simulations. Moreover, we assume that the STNOs are biased with I_dc = 1.5 mA, leading to $ωj(p)/2π≈3.2$ GHz with a coupling strength of $Ωj,j′=25$ MHz (if nodes are connected) and $ωe/2π=2×3.2$ GHz for SHIL.

It has been shown that employing SHIL for annealing in LC oscillator-based IM can increase the probability of achieving the optimal solution.⁹ STNOs exhibit qualitatively different characteristics compared to conventional electronic oscillators, which can be seen from Eq. (3) where the phase $ϕ$ and power p are strongly coupled. Moreover, SHIL in STNOs strongly influences p, which consequently affects the coupling in the network through $pj′/pjΩj,j′$⁠. To explore the effect of SHIL on the STNO-based IM, we employ the annealing schedule shown in Fig. 2(b).

It is worth mentioning that K_e is proportional to the amplitude of the injected RF current at ω_e and consequently, implementing such an annealing schedule in experiments should be relatively straightforward. This annealing schedule was chosen arbitrarily and does not guarantee optimal performance, but comparison of annealing schedules or finding an optimal annealing schedule is out of the scope of this study.

The phase difference between STNOs employing the previously presented annealing schedule is shown in Fig. 2(c). During the first 100 ns, the STNOs distribute relatively uniformly from -π to π as a consequence of the phase coupling $βj,j′$⁠. The state they reach within these 100 ns is already close to the optimal solution, provided that one reassigns their individual continuous phase values to the closest $0°/180°$ phase state. At t = 100 ns, K_e is ramped from 0 to 25 MHz over 10 ns to phase binarize the network, which forces the STNOs to settle to a phase close to $0°$ or $180°$⁠. This process is then repeatedly performed to search for the global optimal/Max-Cut solution. While this small network finds its solution already during the first such step, larger networks require many annealing steps to find their minimum. Every time the SHIL is switched off, there is a substantial delay before the STNOs return to their intrinsic phase state. This delay is longer for stronger SHIL and is a consequence of the sum term in Eq. (3b), which is very close to zero in the SHIL state; the stronger the SHIL, the closer to zero those terms will be.

We then explored the impact of different sizes of Möbius ladder graphs using the annealing schedule in Fig. 2(b). In Fig. 3(a), the probability of achieving the optimal solution in 100 simulation runs with and without SHIL as a function of Möbius ladder graph size is presented. The phases $ϕ$ of all STNOs are initialized with a random uniform distribution between [0, 2π), while the initial power p is kept as 0.3, which is close to the free running STNO power. In the case without SHIL, the final phase difference is rounded to the nearest multiple of π to extract the cut value. Figure 3(a) confirms that SHIL improves the Max-Cut probability in the proposed STNO-based IM. Additionally, in Fig. 3(b), a comparison between the achieved cut as a function of time with and without annealing for a single run for N = 288 is shown. For both cases, the network quickly reaches a solution close to the Max-Cut (within the first 100 ns), but SHIL perturbs the STNO array, leading to a Max-Cut solution, while without SHIL, the system is stuck, close to the Max-Cut. To further confirm this conclusion, the inset in Fig. 3(b) shows the probability distribution of cuts for N = 288 over 100 runs. Without SHIL, solutions close to the Max-Cut exhibit the highest probability, while SHIL modifies the distribution, increasing the Max-Cut probability.

We have verified the functionality of the STNO-based IM and shown that SHIL significantly increases the solution quality. However, Möbius ladder graphs are classified as “easy” Max-Cut problems.⁴⁰ To further demonstrate the potential of the proposed architecture, we employ random cubic (3-regular) graphs which can be considered as “hard” Max-Cut problems.⁴⁰ An example graph of size N = 16 is shown in Fig. 4(a). We generate 20 random graphs of each size, ranging from 16 to 256, and again employ the annealing schedule shown in Fig. 2(b). The Max-Cut solutions are determined using LocalSolver, a commercially available optimization tool, and compared with the solution found by the STNO IM. The results are presented in Fig. 4(b) where each data point represents the average probability of 20 random graphs in 100 runs (total of 2000 simulations). The error bars represent one standard deviation of the probability, which indicates the distribution of hardness for different graphs of the same size. As expected, the probability of finding the Max-Cut solution for random graphs is significantly lower than that for the Möbius graphs, confirming the hardness of these problems. Fig. 4(b) also presents the probability of finding a sub-optimal solution, which is less than 5 cuts from the Max-Cut solution.

It is noteworthy that in all our simulations, we never observed any chaotic behavior, which has been shown to arise in coupled STNOs under certain conditions, in particular, as the number of STNOs increases.^39,41 It was also shown that it may even be beneficial for STNO-based reservoir computing to operate at the edge of chaos.⁴² While it is beyond the scope of our work to address the impact of possible chaos, which could possibly appear in much larger STNO networks, it is something that future studies should investigate.

We have shown that STNOs exhibit fundamental characteristics required to develop an ultrafast and nanoscale IM. However, to achieve an experimental demonstration of a spintronic oscillator-based IM, further research on synchronization of large STNO/SHNO networks is required. Specifically, a general purpose IM requires an all-to-all coupling, which is impossible to achieve using, e.g., dipolar coupling alone. One potential solution might be to employ dipolar/exchange or spinwave coupling for physically close STNOs, while electrical coupling could be used for longer distances. Additionally, an efficient method to control the coupling strength and even the coupling phase needs to be developed to demonstrate a general purpose IM. From a fundamental perspective, the impact of phase noise, variability, and graph density needs to be further explored to fully understand the potential of the proposed spintronic oscillator IM.

This work was supported by the Swedish Research Council (Dnr. 2016-05980) and the Horizon 2020 research and innovation programme (ERC Advanced Grant No. 835068 “TOPSPIN”).