Ising computation based combinatorial optimization using spin-Hall effect (SHE) induced stochastic magnetization reversal Free

Efficient computing models for combinatorial optimization problems have attracted considerable research interest. This trend is in sync with the plenty of information in the present Internet-of-Things (IoT) era. Such huge amount of collected data from multiple sensors is required to be handled properly for certain purposes. As an example, real time management of a smart building (including lighting, cooling, and heating) requires complex optimization using data from multiple sensors.¹ However, solving optimization problems in an efficient way based on conventional computing model is very challenging. Typically, to find an optimum solution for such problems, a performance index has to be computed and compared for every possible input combinations.² However, the computational cost associated with a combinatorial optimization problem increases exponentially with the number of variables. Moreover, if we consider the process of problem solving based on conventional computing (by following a sequential fetch, decode, and execute cycles), finding an optimum (even near optimum) solution seems infeasible keeping in view the energy and power requirements.

Ising model has been researched extensively owing to its simple architecture and inherent ability to solve combinatorial optimization problems.^3,4 A CMOS based implementation of the Ising model can be found in Ref. 5. CMOS implementation, however, requires complex circuits for implementing some of the primitives required for the Ising model like random number generators, majority logic, etc.

As opposed to CMOS implementations, use of spin devices for Ising model opens up new avenues for non-volatile hardware realization of the Ising model. In this manuscript, we propose a non-volatile Ising cell based on controlled stochastic switching dynamics of the magnetization direction in a nanomagnet with an underlying Heavy Metal (HM) layer. The device used in the present manuscript has been demonstrated experimentally in Ref. 6; moreover, a simulation study presenting a random number generator and an artificial neuron can be found in Refs. 7 and 8, respectively.

Recently, an Ising model based on a similar spin Hall effect (SHE) based device and low Energy Barrier (EB) (<5 KT) was proposed in Ref. 9. However, the use of magnets with low EB severely limits the design space, since it has to be ensured that the read current is small enough to avoid possible pinning of the free-layer during the read operation. The constraint on the read current and hence the read voltage also make it challenging to drive multiple devices from a given read port of a particular nano-magnet. Further, the use of low EB magnets compromises the intrinsic non-volatility of spin devices. Another proposal, similar to Ref. 9, but based on high EB magnets initially biased towards the magnetic hard axis can be found in Ref. 10. The inclusion of dipole coupled magnets in the proposed device presented in Refs. 9 and 10 increases the complexity of the device. Besides, details about circuit requirements for interconnections and weighted summation of currents are missing.

The simplicity of the device used in the present work in addition to the presented CMOS circuits allows us to construct an Ising cell, which can be seamlessly interconnected to solve combinatorial optimization problems. A numerical simulation framework based on the stochastic Landau-Lifshitz-Gilbert (LLG) equation is developed to analyze the switching characteristic of the proposed device. Further, by solving coupled stochastic LLG equations and SPICE simulations, we show solutions for some representative combinatorial problems obtained by using our proposed Ising cell.

Before we describe the proposed device, we would give a brief introduction of the Ising spin model, in general. The Ising model considers the behavior of magnetic spins and the coupling between them. Figure 1(a) illustrates a simple view of the Ising model and the definition of associated Hamiltonian (H)—the total energy of the system. The model consists of individual spin state (s_i), interconnection coefficient between two spins (J_ij), and external magnetic field (h_i). Each spin can have one of the two states, up and down, and there are four interconnections with neighbors in this model. The spin at the center is named as s_C and the four neighbors are s_U, s_D, s_L, and s_R. The interconnection weights between s_C and its neighbors are denoted as J_CU, J_CD, J_CL, and J_CR, respectively. These weights model the coupling between spins and are used to determine next state of the spins. For example, if J_CU has positive sign (i.e., +1), it implies s_U tries to align s_C parallel to itself. Likewise, a neighboring spin with a negative weight tries to align the given spin anti-parallel to itself. Since there are four neighbors, the next state of s_C is decided based on a majority vote—if majority of the neighbors of a given spin state want to keep the given spin in +1 direction, then the next state of s_C will be +1, else it would be −1. The Hamiltonian (H) also changes as the states of the spins are updated. Hence, once the problem is mapped to the system properly (by programming the weights for each interconnection), the system tries to reach the energy minimum state by switching the states of the spins through the aforementioned majority coupling. When the system reaches the global minimum energy state, the solution is obtained by examining the final states of the spins.

Figure 1(b) shows the total energy (H) of the system as a function of different spin states. As shown in the figure, the energy profile has a global minimum and also multiple local minimum states. This implies that the system could easily get stuck at the local minimum state during the process of problem solving (i.e., the system evolves to different states through coupling). To avoid the system being stuck into a local minima, annealing process, such as “Simulated Annealing (SA)”^11,12 and “Quantum Annealing (QA),”^13–15 has been proposed. During a SA process in conventional Ising model, the system starts from a known initial state at a non-zero temperature. The system then evolves towards the minimum state of the Hamiltonian by lowering its temperature gradually. In contrast, in a QA, the temperature can be replaced by a quantum mechanical effect such as probabilistic quantum tunneling.¹⁴ Whether it is QA or SA, annealing always includes some kind of randomization of the next state logic to get the system out of the local minima.

Despite the fact that SA and QA can find the lowest energy state of the Ising spin model efficiently, the implementation of such a system needs control of the state of each spin and coupling between them. Also, the state of individual spins needs to be monitored for total energy of the system which is challenging from hardware perspective. This is why hardware implementation of the Ising model did not receive much attention, even though theoretical background has been widely explored by the research community. Recently, hardware implementations of the Ising spin model have been proposed using superconducting material¹⁶ and CMOS only implementation.⁵ In CMOS circuits, the annealing process can be implemented by generating a random address that is used to choose a specific spin to be flipped.⁵ However, such implementations require complex hardware for random number generation, address decoding, and write operation for the specific spin state. These series of operations have to be executed multiple times to get the system state out of the local minima. Furthermore, randomizing spin states based on random number generation can potentially move the system to totally different state so that the system might not reach the global minimum state eventually.

Another difficulty in the implementation of the Ising model is due to the complexity of the majority voting circuit. As explained, the next state of each spin is determined by the interactions with all neighboring states. Multiple solutions might be possible to implement majority function based on digital, analog, and even with mixed-mode design. However, all of these implementations are expensive in terms of silicon area and power consumption due to its multi-input nature and complexity of operation.

The motivation of our research arises from the fact that a spintronic device, such as a Heavy-Metal (HM) based Magnetic-Tunnel-Junction (MTJ), can potentially provide the aforementioned two important characteristics required for the Ising model, viz., (1) stochasticity (random spin flip required for annealing) and (2) majority voting. In order to mimic the aforementioned functionality, a nanomagnet is required that is able to switch its states with a certain probability to facilitate the annealing process and also change its state based on a majority vote which is crucial for the system to evolve towards minimum energy state. In Sec. II, we describe the mapping of the magnetization dynamics of a nanomagnet to the Ising operations such as a natural randomizer and a majority vote logic.

Let us first describe the basic device structure used for our Ising model and its principle of operation. Figure 2 shows a three terminal device structure, consisting of the conventional MTJ stack formed by a Tunneling Barrier (TB) sandwiched between two nanomagnets. Since the magnetization of the upper ferromagnetic layer is fixed, it is termed as the Pinned Layer (PL). On the other hand, the magnetization of the bottom layer, denoted as the Free Layer (FL), can be manipulated by an incoming spin current. Depending on the direction of the FL, the MTJ structure can have two stable states. If the magnetizations of the two ferromagnetic layers are in the same direction, it is in the parallel configuration (P); otherwise, it is in the anti-parallel (AP) configuration. These two states have different resistances across the vertical direction of the device. Typically, the AP state has a higher resistance, and the ratio between P and AP states is defined as the Tunnel MagnetoResistance (TMR) ratio.

The FL is in contact with a heavy metal (HM) like Pt or Ta. On passing a charge current through the HM, from terminal T1 to T2, the direction of magnetization in the FL becomes parallel to that of the PL. When the direction of charge current flowing through the HM is reversed (from T2 to T1), the FL switches its direction such that it is now anti-parallel to the PL. The switching of FL due to a charge current flowing through the HM can be attributed to the large spin-orbit-coupling (SOC) of the HM. SOC is a relativistic effect that theoretically originates from the full relativistic wave-equation. In the current scenario, due to SOC based effects, like the Spin Hall Effect (SHE), the electrons flowing through the HM are deflected such that up and down spins split. As shown in Fig. 2, the up-spins are deflected in the +z direction and down-spins in the –z direction. This spin splitting results in a resultant spin current flowing in the +z or –z direction, based on the direction of the charge current. The spin current (I_s), thus generated due to the charge current flowing in the HM, exerts a torque on the magnetization direction of the FL, making it parallel or anti-parallel to the PL. Thus, a charge current flowing between the terminals T1 and T2 can switch the state of the MTJ from R_AP to R_P or vice-versa. In order to sense the resistance of the MTJ, a voltage can be applied between terminals T3 and T1/T2. By sensing the current flowing through the MTJ stack (or detecting the voltage level across the device), one can conclude if the current state of the MTJ is R_AP or R_P. Based on the above description, the three terminal structure of the device, as shown in Fig. 2, offers the following desirable characteristics: (1) the write and read path are decoupled and can be independently optimized; (2) the efficiency of spin current generated from the charge current flowing through the HM can be greater than 100%,¹⁷ resulting in low write-energy; and (3) by controlling the amount of current flowing through the HM, one can not only alter the switching probability but also implement a majority function. Later, in the manuscript, we would describe how these desirable characteristics of the HM based MTJ device can be used efficiently to mimic the various operations required for the Ising spin model.

We would now describe the equations used for modeling the device as shown in Fig. 2. Under an external excitation, for example, a magnetic field or a spin current, the magnetization dynamics of the FL can be obtained by the well-known Landau-Lifshitz-Gilbert (LLG) equation under the single domain approximation with an additional term for the spin transfer torque proposed by Slonczewski¹⁸ 

where $m̂$ is the unit vector corresponding to the direction of magnetization in the FL, $γ=2μBμ0/ℏ$ is the gyromagnetic ratio for electron, α is Gilbert damping ratio, H_eff is the effective magnetic field including the shape anisotropy field¹⁹ and uniaxial interface anisotropy field,²⁰ $Ns=MsV/μB$ is the number of spins in the FL volume V (M_s is saturation magnetization and μ_B is Bohr magneton), I_s in Eq. (1) is the spin current flowing through the FL in the +z or –z direction, and $M̂$ is magnetization direction of PL.

Based on the recent experimental studies,^6,21–24 the spin current generated due to a charge current flowing through the HM can be estimated as $Is=θSHWMTJtHMIQ$⁠, where I_Q is the charge current flowing through the HM, θ_SH is the spin-hall angle,²² W_MTJ and t_HM are device dimension parameters, as shown in Fig. 2. The details of the device parameters we used for benchmarking are summarized in Table I.

Finally, to model the effect of thermal noise at non-zero temperature, the thermal noise is accounted as a random thermal field,²⁶ $Hthermal=α1+α22KBTγμ0MsVδtG0,1$⁠, where $G0,1$ is a Gaussian distribution (zero mean, unit standard deviation), K_B is the Boltzmann constant, T is the temperature, and δ_t is the time step. Under the effect of thermal noise, the switching behavior of the MTJ during the write operation due to the charge current flowing through the HM layer becomes stochastic in nature. Also, the probability of switching changes according to the magnitude of the input charge current. Figure 3(a) illustrates a graph showing switching probability (P_SW) versus input charge current (I_q) through the HM layer. The charge current pulse was applied for a fixed time (3 ns), and an additional 6 ns wait time was included for the magnetization to relax. The amount of charge current varies from 40 μA to 160 μA with 5 μA step, and 10⁴ simulations were executed for each simulation step to get a switching probability. As shown in the figure, when an input current of 90 μA is applied for 3 ns, the P_SW becomes roughly 50%. Figure 3(b) shows the normalized magnetic moment change with 100 write operations (assuming that the initial magnetization is 1 and is being flipped to −1 direction) when the input charge current is kept at 90 μA for 3 ns followed by 6 ns of relaxation time. It can be observed that 51 out of 100 cases flipped, which is close to the ratio one would expect from Fig. 3(a).

The stochastic switching of the nano-magnet has been attracting significant interest in recent years and is being used for intriguing applications such as neuromorphic computing.^27,28 In this paper, we use this stochastic flipping nature of the nanomagnet as a natural randomizer—one of the key elements for the “natural annealing” process. The baseline idea of a general annealing process in Ising model lies in perturbing the spin states randomly to get the system out of the local minimum energy state. Thus, while switching the state of the MTJ, we can tune the input charge current flowing through the HM such that the MTJ switches with the desired probability. The write current can be controlled with ease by adopting simple CMOS peripherals which will be explained later. It is worth noting here that the proposed natural annealing can also provide time-varying switching probability (by adjusting the input current value) which mimics the natural property of annealing through the temperature control. This helps to find a global minimum quickly when the system reaches the end of iterations.

Additionally, we can exploit another benefit from the nanomagnet due to the strong dependence of the switching process on the input charge current flowing through the HM. As explained, charge current flowing through the HM layer induces spin current in transverse direction at the FL-HM interface, which flips magnetization of the FL. The dependence of the flipping process on the input current can be explained by the energy profile of the FL. Assume that the angle between FL magnetization and the PL magnetization is represented by θ. The FL energy as a function of θ is shown in Fig. 4(a), where the two stable states (θ = 0° and θ = 180°) are separated by the energy barrier E_B. Here, we assume that the MTJ changes its state from P (θ = 0°) to AP (θ = 180°) state. Once the input charge current is presented to the HM layer for a duration t_WRITE, spin current is induced based on the spin-hall effect (SHE) and makes the spin at point 1 to move uphill along the energy profile. If the charge current is not enough to move the spin across the energy barrier, the spin stops at the metastable state (point 2) and falls down again to the point 1 during t_RELAX. On the other hand, once enough charge current is applied, ultimately the spin will overcome the energy barrier and move to the point 3. Figure 4(b) shows these situations using magnetic moment change and two current pulses with different magnitudes. Note that, in the presence of thermal noise, a hard switching threshold current does not exist. A particular amount of current can only flip a magnet with certain probability unless the applied current is large enough to deterministically switch the magnet.

The decoupled write operation in HM based MTJ can be used as an efficient way to construct the majority vote function required for the Ising model. Let us consider a given connection between a particular spin state and one of its neighbors. We would represent the connection as (+1,+1), where the first number denotes the spin state (+1 represents up spin and −1 down spin), while the second number in the bracket represents the associated weight of the given spin and its neighboring spin. As mentioned in Sec. I, the next state of each spin is determined through the interaction with all connected neighboring spins. For example, if a neighbor has a state (+1,+1) or (−1, −1), then it would want the next state of the spin to be +1 (obtained by product of spin state and associated weight). In case of (+1, −1) or (−1, +1), the neighbors would tend to make the next state of the given spin under consideration as −1. In the presence of multiple neighbors, a majority vote is taken to determine the next state of the spin under consideration. The hardware implementation of such a multiplication (product of spin state and associated weight) and majority vote functionality requires complex circuit.

Our proposed HM based MTJ circuit that can efficiently implement the majority vote functionality is shown in Fig. 5(a). The HM layer receives the input current proportional to the number of voters from the neighboring spin states using multiple current sources and switches. The corresponding P_SW is depicted on the top of Fig. 5(b). Here, we assume the possible current range for write operation is limited within 60 μA to 120 μA. This range is equally divided among its neighbors. For instance, if there are four neighbors, each neighbor would contribute 15 μA of write current by voting to one of the two potential states. Note that the leftmost current source is used to provide a default magnitude of current (60 μA) during the write cycle. Based on this, if there are 0 voters, the current during the write operation becomes 60 μA, which leads to 2% P_SW. If there are 4 voters, then the current becomes 120 μA and corresponding P_SW becomes 96%. This directly mimics the general rule of majority vote—low P_SW with less voters, high P_SW with more voters. Note that, due to its probabilistic nature, there are chances of spin flip into unwanted state. This is not an issue, it rather mimics the natural annealing process as discussed earlier.

The overall device-circuit configuration for single spin model is shown in Fig. 6. For read operation, reference resistor (R_REF) and the switch transistor are used in series on the top of the nanomagnet. The resistance of R_REF is in between two stable resistance states of the MTJ device. Hence, if the MTJ resistance is smaller than R_REF (in parallel magnetization configuration), the output of the inverter becomes high and vice versa. The output of inverter, i.e., the next state of the current spin is then stored in the following D-latch. This additional storage element de-couples the evaluation of the next state from interacting with other neighboring cells. Consequently, every Ising cell can be updated concurrently at each operation step regardless of the size of the system, thereby improving the scalability of the proposed architecture. The “next state” signal is used to choose the direction of the current flowing during the write operation and also be presented to neighboring cells for the majority vote operation. For the majority vote and write operation, there are current sources with multiple branches (each with stacked transistors, one for biasing and the other for switch operation) and XNOR gates to combine information from the neighbors (product of spin state and interconnection weights). Note, here we assumed that there are four neighbors. The number of neighbors can be extended by adding more branches on current source and controlling the amount of current from each branch. Thus, (1) the XNOR gate (which implements the multiplication function), (2) the transistor switches and the dependence of P_sw on input current (which mimics the majority vote logic), (3) the “next state” logic which controls the direction of current flow, and (4) the inherent stochasticity of switching (representing the annealing process) constitute our basic Ising cell, which can be replicated to implement combinatorial optimization problems.

In this section, we would first elaborate how the basic Ising cell described in Sec. II, can be used to map combinatorial problems. For the Ising spin model shown in Fig. 7(a), each spin state consists of a MTJ device with underlying HM layer and associated interface circuits. The interconnection weights between neighboring spins are a function of the specific problem to be solved. These weights are programmed into the system before the Ising model tries to converge towards the minimum energy state. Based on the initial spin states and associated weights, the system evolves by updating the state of each spin through the coupling (annealing and majority vote) described earlier in the manuscript. The time evolution of the Ising system in our simulation framework was achieved by following the steps as shown in Fig. 7(b). After selecting a specific spin to be updated, the neighboring states and the weights associated with each interconnection are examined. Then, the number of voters for the potential next state (either +1 or −1) and corresponding charge current for write operation is obtained through SPICE simulations. The resultant current level is then fed to the LLG solver to analyze the magnetization switching behavior of the HM based MTJ. After applying the current pulse for a duration of 3 ns to the HM layer followed by 6 ns of relaxation time, the final direction of the FL magnetization is examined, which determines the next state of the spin under consideration. This process is repeated until the state of the all the spins in the system is updated. Then, the system computes the Hamiltonian of the next state to check whether the system has found a solution.

The aforementioned generic methodology was the applied to two practical combinatorial optimization problems—Maximum-cut and Graph coloring. First, the proposed Ising spin model is used to solve Maximum-cut problem. The problem can be defined as finding two mutually exclusive subsets of spins by connecting edges (which connect spins from two separate subsets) so as to maximize the summation of weights along the edges (i.e., boundary between two subsets).²⁹ We have used ∼3900 spins for this application and the interconnection weights are programmed so that the spin states show the digit numbers from 0 to 4 without noise once the system reaches the minimum possible energy state. Figure 8 shows the results of the maximum-cut problem and also transition of the system energy over iterations. As shown in Fig. 8(a), the system energy shows a rapid decrease up to ∼100 iterations followed by a slow decrease in energy for additional ∼300 iterations. The initial rapid energy drop happens due to the changes in the spin states from the initial random states to a set of two separate regions with some associated noise. This implies that the system reaches local minima. Then the system seeks for the solution with global minima based on the natural annealing. Note that this global minima search process could be expedited through the time-varying switching probability function obtained by controlling the current flowing through the HM underlayer. Here, the amount of current presented to each Ising cell was adjusted at 400th iteration so that the system have less random flipping, which leads to faster convergence. This directly mimics the conventional SA process, wherein lowering the system temperature is represented by the decrease in the switching activity of the spin states. Figure 8(b) shows visualized spin states at the initial, intermediate, and final steps. Here, the black dot denotes spin state −1 and white dot represents spin state +1. Initially, spin states start from random distribution of −1 and +1 (black and white dots). As the states are updated through the coupling with neighboring states, target digit numbers are visible with some noise (100th iteration). After 400 iterations, the system almost finds the solution, but, still there exists nontrivial amount of noise mainly due to the “natural annealing.” The clear output image shows that the optimum solution of the current maximum-cut problem has been obtained.

Next, our proposed Ising spin model with HM based MTJ was applied to the Graph coloring problem. Graph coloring is a famous NP-complete problem³⁰ and is defined as “Is it possible to color n-vertices with k-colors such that no two adjacent vertices have the same color?” To implement specific hardware to solve the graph coloring problem using our Ising spin model, a pre-processing step is needed to prepare a weight matrix. This can be generated according to the penalty Hamiltonian in Ref. 31. This Hamiltonian basically defines rules to be obeyed and applies an energy penalty whenever these rules are violated. For example, one term of the Hamiltonian for this particular problem provides energy penalty once the edge connects two vertices with the same color. Consequently, optimum solution can be obtained when there is no energy penalty, hence the system evolves towards minimum energy state. Interested readers are directed to Ref. 31 for more details on the Hamiltonian for the Ising model. Once the weight matrix is prepared, it shows interconnect map and also weights for each interconnections. We prepared 3 simple Graph coloring problems and implemented corresponding hardware based on weight matrix from penalty Hamiltonian. Figure 9 shows the details of the problem from our proposed Ising solver. It is worth noting here that since the spin can only have one of the two states, a total of n × k spins are needed to represent all possible states for the graph coloring problem (with n-vertices and k-colors). In this case, each spin state is denoted as v_i,j, where i represents current vertex and j represents current color. For example, v_1,1 can be interpreted as a spin representing vertex 1 and color 1. The simulations were conducted 1000 times for each problem to get an average number of iterations to reach minimum energy state. The transition of the system energy is monitored by checking the states of all the spins after each epoch to check if the system has found a solution.

Finally, let us briefly discuss the energy consumption of the HM based MTJ device used in the Ising spin model. The operation of a single Ising cell can be divided into three parts—read operation, write operation, and relaxation time. The time duration for write, relax, and read cycle was taken to be 3, 6, and 1 ns, respectively. The energy consumption during the write cycle is mainly due to the input charge current flowing through the HM layer. Assuming an average input current of ∼90 μA, the energy consumption during the write operation is evaluated to be ∼0.27 pJ with VDD of 1 V. Our simulations indicate that the dynamic energy consumption per cell due to CMOS peripherals is small as compared with the energy dissipated in the MTJ during the write operation. Likewise, the device-circuit simulation of the read circuit including the associated CMOS transistors yielded an average energy consumption of ∼0.04 pJ, when considering the average read current to be ∼38 μA and VDD to be 1 V. In contrast, the energy consumption involved in the relax mode and CMOS switches resulted in insignificant contribution (∼0.01 pJ) to the total energy consumption. Overall, the proposed single spin model based on our HM based MTJ along with peripheral circuits consumes ∼0.32 pJ of energy per single spin update operation. Even though the write operation consumes most of the energy required (∼83%), we believe with improvement in material parameters, write energy can be significantly reduced.

In summary, we have proposed SHE-MTJ based Ising cell to solve combinatorial optimization problems. We demonstrate the mapping of annealing and majority vote functions to the behavior of the spins in the nanomagnet. Although the stochastic switching nature of the MTJ due to the thermal noise is usually regarded as a problem in typical memory applications,³² such random switching can be exploited to build Ising spin model having the property of “natural annealing.” Also, the decoupled write and read path through the HM underlayer enables the majority vote function with minimum number of devices. Using coupled magnetization dynamics and SPICE simulations, we present solutions for two combinatorial optimization problems—Maximum-cut and Graph coloring. We believe that the behavior of our HM based MTJ device, mimicking the key elements of the Ising spin model, can potentially pave the way for scalable, low-power, and simple Ising solver capable of handling complex combinatorial optimization problems.

The work was supported in part by Center for Spintronic Materials, Interfaces, and Novel Architectures (C-SPIN), a MARCO and DARPA sponsored StarNet center, the Semiconductor Research Corporation, the National Science Foundation, Intel Corporation, and the National Security Science and Engineering Faculty Fellowship.