Magnetoelectric oxide based stochastic spin device towards solving combinatorial optimization problems

Solving combinatorial optimization problems is challenging. Mapping onto the ground-state search problem of the Ising Hamiltonian is a promising approach in this field, where the components of the optimization set are modeled as artificial spin units. The search for a suitable physical system to realize these spin units is an active area of research. In this work, we have demonstrated a scheme to model the Ising Hamiltonian with multiferroic oxide/nanomagnet units. Although nanomagnet-based implementation has been shown before, we have utilized the magnetoelectric effect of the multiferroics to make voltagecontrolled spin units with less current flow in the network. Moreover, we have proposed a unique approach of configuring the coupling network of the system directly from the Ising Hamiltonian of a traveling salesman problem (TSP). We have developed a coupled micromagnetic simulation framework and solved TSPs of size 26-city and 15-city with an accuracy of 100% for the latter.

salesman problem (TSP), leading to huge current flow through the whole network. For example, if N = 10, the total current requirement roughly amounts to 100 × I c , where I c is the current through individual spin unit. This poses challenge to the scalability of the system as the problem size increases. On the other hand, in our device, the write unit consists of a multiferroic oxide/nanomagnet stack, where the magnetization state is controlled by the voltage applied across the oxide. This oxide is thick enough to inhibit current flow in all of the N 2 write units.
• Read: Reading the magnetization state in both methods are performed by a current flow through an MTJ, where the free layer is coupled to the nanomagnet of the write unit through dipolar interaction. However, the output parameter in previous works is current (proportional to the resistance of the MTJ), which is passed on to the inputs of interacting spin units. But, in our device, the MTJ current is converted to an output voltage by using a voltage divider circuit. This voltage later provides inputs to coupled units. • Interconnect: The spin-spin interaction is emulated by sending the output from one spin unit to the other.
The current-based input-outputs in previous works cause current flow from one interacting unit to the other, whereas, in our device, there is no current flow between connected spin units due to voltage-based coupling. This is an important aspect of scalability of the scheme. Also previous works suggest additional amplification 9 or CMOS circuits 10 in between coupled spin units in order to drive the fan-out, whereas in our device, direct cascading is possible by proper device engineering.
It is to be noted that the network of spins in our model does not entirely resemble a Boltzmann machine 11 . The focus of our scheme is to make the lower energy states to be more probable with no restriction on the overall energy landscape. Therefore, although the probability of observing a state with energy E gradually increases with decreasing E, it may or may not follow the Boltzmann distribution.
In this work, we have presented the structure and switching characteristics of these voltage-controlled spin units. Additionally, we have showed a way of direct configuration of the spin-spin and spin-external source coupling network from the Ising Hamiltonian, taking TSP as an example. As a validation of our model, we have developed a coupled LLG (Landau-Lifshitz-Gilbert) equation-based simulation framework and solved TSP problems of size 15-city ((N − 1) 2 = 196 spin units) and 26-city (625 spin units). We compared our results with existing heuristic (Lin-Kernighan) algorithms 12 used in the field of computer science and achieved 100% accuracy for 15-city.

Ising formulation of Traveling Salesman Problem
Ising model, although originally introduced as a model for ferromagnetic materials, is widely applied in molecular biology, chemistry and other areas due to its combinatorial interpretation. According to this model, the spin dynamics in a ferromagnetic lattice, consisting of N lattice sites, is governed by the following Hamiltonian 2 : Here x i is the spin state of the molecule at the i-th lattice site, which can assume either 'up' or 'down' state. J i,j and h i correspond to the energies due to the interactions with the nearest neighbors and external fields, respectively. The spins interact in such a way that they tend to eventually line up in the configuration producing the lowest value of H, thus transitioning from a high energy random state to a low energy ordered state at or below a critical temperature. Finding the spin configuration which minimizes H is itself an NP-hard problem. Hence, the elements of a combinatorial optimization problem can be thought of as a collection of spins, x i , where the Ising energy function, H represents the parameter to be optimized. For example, the traveling salesman problem, which asks for the ordering of cities to visit so that the total distance travelled is the minimum, can be formulated as an Ising energy function in the following way 13 : For an N city TSP, the system consists of N 2 bits/spins x u,i ( Fig. 1(a)), which can take either '0' (down) or '1' (up) depending on the fact whether city u will be travelled at order i or not. W uv is the distance between cities u and v.

Nanomagnet-based implementation
The fundamental building blocks of the Ising model, also known as artificial spin units, act as a random number generator (RNG) (randomly switches between 'up' and 'down' states) unless acted upon by external force or neighboring interactions. The stochastic switching characteristic of nanomagnets, arising from the inherent thermal noise, enables its use as RNGs. Nanomagnet-based RNGs have been demonstrated previously, like spin transfer torque based spin-dice 14 , spin-orbit torque based spin-dice 15 , voltage controlled spin-dice 16 etc. However, high current requirement in spin-torque-based devices and unipolar switching 17 of VCMA-based devices have stirred the quest for new mechanisms, like the magnetoelectric (ME) effect.
Device description. Among different ferromagnet/ME oxide heterostructures 18 , BiFeO 3 (BFO)/CoFeB is found to demonstrate exchange bias coupling, which is suitable to accomplish 180° switching 19,20 . Exchange bias interaction occurs at the interface between the ME oxide and the ferromagnet, while the spin and charge polarization of the multiferroic material BFO are coupled to each other. Hence, by switching the electric polarization of BFO, the magnetization of CoFeB can be switched. In our model, we have used an ultrathin CoFeB film (thickness of 0.9 nm) with perpendicular magnetic anisotropy (PMA) as the ferromagnet in contact with a multiferroic layer BiFeO 3 (BFO) which consists the write unit ( Fig. 1(b)). Each spin x u,i is represented by the magnetization moment of CoFeB. The application of an electric field, V IN across the ME oxide creates an effective magnetic field, → B ME experienced by x u,i . Magnetization reversal of x u,i takes place depending on the strength and direction of B ME → . Electrical switching of this type of PMA CoFeB/BFO system has been experimentally demonstrated in ref. 20. However, since the crucial aspects of these fields, like coupling mechanisms, switching dynamics etc. are still under intense research, instead of following any particular experiment set, we have used a generic parameter called magnetoelectric coefficient 21 (α ME ) to model the ME effect in our device. (The details of the simulation methodology is explained in the next section.) The thickness of the BFO layer is 5 nm to inhibit tunneling current in the write unit. The area of the CoFeB layer used here is small enough (16 nm × 8 nm) to make it work as an RNG when unbiased. Figure 1 The read unit ( Fig. 1(b)) contains a voltage divider circuit consisting of an MTJ and a reference resistance, R REF . This read circuit is electrically isolated from the CoFeB/BFO stack by an insulating layer. However, dipolar interaction couples the free layer of the MTJ to the magnetization x u,i of the CoFeB layer in the write unit. The resistance of the MTJ, R MTJ varies between parallel (R P ) and antiparallel (R AP ) configuration depending on the value of x u,i and the output voltage, V OUT changes accordingly.
Here, V REF denotes the supply voltage. The variation of V OUT with the magnetization moment is demonstrated in Fig. 2(b), generated by a behavioral model (explained in the next section). It is to be noted that V REF should not exceed the switching threshold of the MTJ free layer to keep its state unperturbed. At the same time, it should be sufficient to let V OUT reach the desired values to drive subsequent units. Hence, optimization of the MTJ dimensions (i,e, R MTJ ), R REF and V REF is necessary. It is worth mentioning here that the electrical isolation of the READ and WRITE unit demands two free layers for each section. It helps prevent any undesired perturbation of the current magnetization state of the WRITE free layer by the READ current.

Neighboring and external interactions. The input voltage V IN consists of two components, V neighbor
and V external which model the interaction with neighboring units and external sources, respectively. In other H A imposes the restriction that each city should be visited only once. It can be reformulated into the following form: A v j u i , ,

Each term in H
Each term in H B dictates whether city u and v should be visited consecutively depending on the distance between them, W u,v . At the ground state of the system, H B refers to the minimum distance travelled. With this restriction in mind, the truth table in Fig. 3(a) works towards selecting the spins in consecutive columns based on their distances. Note that, in our model, a spin x v,j in state '1' implies that city v is chosen to visit at order j. Each spin x v,j from Equation 4 is situated at (i + 1)-th or (i − 1)-th columns, whereas x u,i lies at the i-th column ( Fig. 3(a)), given the columns dictate the order of travel. Each term in Equation 4 is modeled as an external input voltage V u,v being applied to the spin unit x u,i , while this voltage source is switched 'ON'/'OFF' by the states of the spin units x v,j in the adjacent columns ( Fig. 3(b)). In case of the last column (column N), the adjacent columns are (N − 1)th  (Fig. 3(d)). Hence, we are trying to determine which city among B,C and D should be visited at order 2, given city A is visited first (x A,1 = 1). The switches are always 'ON' (because x A,1 = '1'). In addition, all units are being injected with V neighbor , as described in the previous paragraph, to prevent simultaneous selection of spins in the same column. In this simulation, we have chosen W A,B < W A,C < W A,D . Therefore, Fig. 3(e) demonstrates that '100' is the favorable state for x B,2 x C,2 x D,2 , implying that city A to B is the favorable path, not A to C or A to D. Also, H B is centered at W AB (Fig. 3(e)), which is the shortest route.

Simulation Methodology
In order to solve an N by N problem, we have used an N − 1 by N − 1 network of spins keeping node 1 fixed to appear first in the cycle. To calculate the magnetization dynamics of these (N − 1) 2 nanomagnets, we have developed a simulation framework consisting of a set of (N − 1) 2 stochastic LLG equations coupled with each other through their input voltages, V IN . In our numerical simulation, we have used a time step of 0.02 ns. After each time step, the magnetization moment of each bit is updated, and these updated values are used to generate the input voltages for the next time cycle (Fig. 4). A behavioral model is used to generate the input voltages, taking corresponding spin values as inputs. The details of the micromagnetic simulation and the behavioral model are described in this section. Here m is the normalized magnetization moment of the ferromagnet layer, α is the Gilbert damping constant, q is the charge of an electron and γ denotes the gyromagnetic ratio. → H eff is the effective magnetic field acting on the magnetization, which consists of the anisotropic field, → H anisotropy , demagnetization field, → H demag , thermal field, → H thermal and effective magnetoelectric field, → H ME . Note that there is no spin transfer torque due to negligible current flow through the oxide. k B is the Boltzmann constant, T is the temperature of the system, M s is the saturation magnetization, Vol is the volume of the free layer and dt is the discrete time step used in the numerical simulation. ξ → is a 3-component vector whose components are zero mean Gaussian random variables with standard deviation of 1.
We have modeled the ME field through the magneto-electric coefficient α ME The experiment in ref. 27 shows a nonlinear variation in α ME versus the applied voltage, with the values ranging from ~0.001×10 −7 sm −1 to 1 × 10 −7 sm −1 . However, for the voltage range in our model, we have used an average value of 0.03 × 10 −7 sm −1 (1/c, c = speed of light). The parameters used in the LLG simulations are listed in Table 1.
Behavioral model. In order to calculate V IN (in Equation 7) applied to each ME oxide/nanomagnet unit, we have used behavioral models to calculate the constituents V neighbor (Fig. 2(b)) and V external (Fig. 3(b)) which require relevant spins, x u,i /x v,j (generated by LLG equations) as inputs. Since R MTJ follows a sigmoid function with its magnetization moment (SPICE model of the MTJ: ref. 29), V neighbor (i.e., V OUT , which is linearly dependent on R MTJ ) is modeled as A tanh (Bx u,i ), where the parameters A and B are adjusted to fit the desired values of V 1 and V 2   Table 1. Summary of the parameters used in the simulation.
( Fig. 2(b)). On the other hand, V external is linearly related to the distance matrix by the equation: Here W is the value of the corresponding weight/distance, W min and W max are the minimum and the maximum value in the distance matrix, respectively. This behavioral model and the set of coupled LLG equations are solved self-consistently with a time step of 0.02 ns.
In addition, we have used simulated annealing in order to gradually move towards the global minima. The total time required for the nanomagnet array to reach a steady state starting from a random distribution is around 10 μs with an annealing period of 50 ns. At i-th period, |V 2,i | = |V 2,initial | + i × δV neighbor , |V low,i | = |V low,initial | + i × δV external and |V high,i | = |V high,initial | + i × δV external . (V 2 is shown in Fig. 2(b), V low and V high are shown in Fig. 3(c)). In this way, we move towards a steady state as time goes by. The efficiency of the annealing schedule largely depends on the initial points and the cooling rate. In our simulation, we have used |δV neighbor | = 50 mV and |δV external | = 100 mV.

Results
We developed a coupled LLG equation-behavioral model simulation framework to solve TSPs of two sizes: 15 by 15 ((N − 1) 2 = 196 nodes) and 26 by 26 (625 nodes). Figure 5(a,b) show a side-by-side comparison of the routes obtained by brute force search (red) and our approach based on ME oxide/nanomagnet stochastic device (blue), where the case for 15-by-15 matches 100%. However, as the problem size increases, the solution becomes more susceptible to the efficiency of the annealing schedule. Hence, the accuracy goes down for the 26 by 26 problem. However, the aim of our work is not to analyze optimized annealing schedule, rather demonstrate that near-optimal solution is achievable with this ME oxide/nanomagnet based stochastic device. The problem dataset and brute-force solutions are taken from refs 28 and 30. The plots of the Hamiltonian function in Fig. 5(c,e) indicate that the system moves towards low energy ordered state over time as we approach solution. The initial and final magnetization states are shown in the insets. Figure 5(d) shows a comparison of our results with LK 12 heuristic algorithm. It is possible to boost up the accuracy of the 26-city problem by adjusting the annealing schedule. The code for solving the LK algorithm with our test data has been taken from ref. 31.
We have also calculated the amount of energy dissipation in the spin units. As stated earlier, current flow in the write unit and the interconnect is negligible. Hence, power dissipation takes place only in the read unit. Note that the area of the MTJ is large (96 nm × 72 nm) enough to make the critical switching current much larger (high energy barrier) than the current flow in the voltage divider circuit. With these parameters, the maximum power dissipated in the individual spin units amount to 0.15 mW. It can be reduced by using better multiferroic oxides with strong magneto-electric effect, i.e. large magnetoelectric coefficient α. To this date, the value of α (at low voltages) obtained experimentally for exchange bias effect is limited to 1/c to 0.01/c, c is the speed of light 19 .