# Toffoli gadget for magnetic-tunnel-junction Boltzmann machines

Dairong Chen, Augustin Couton Wyporek, Pierre Chailloleau, Ahmed Sidi El Valli, Flaviano Morone, Stéphane Mangin, Jonathan Sun, Dries Sels, Andrew Kent

► To cite this version:

Dairong Chen, Augustin Couton Wyporek, Pierre Chailloleau, Ahmed Sidi El Valli, Flaviano Morone, et al.. Toffoli gadget for magnetic-tunnel-junction Boltzmann machines. Physical Review Applied, 2025, 23 (3), pp.034044. 10.1103/PhysRevApplied.23.034044. hal-05274057

HAL Id: hal-05274057

https://hal.science/hal-05274057v1

Submitted on 23 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A Toffoli Gadget for Magnetic Tunnel Junctions Boltzmann Machines

Dairong Chen,¹ Augustin Couton Wyporek,¹,² Pierre Chailoleau,² Ahmed Sidi El Valli,¹

Flaviano Morone,¹ Stephane Mangin,² Jonathan Z. Sun,³ Dries Sels,¹,⁴ and Andrew D. Kent¹

¹Center for Quantum Phenomena, Department of Physics,
New York University, New York, NY 10003 USA

²Université de Lorraine, CNRS, Institut Jean Lamour, F-54000 Nancy, France

³IBM T. J. Watson Research Center, Yorktown Heights, NY 10598, USA

⁴Center for Computational Quantum Physics, Flatiron Institute, New York, NY, USA

Magnetic Tunnel Junctions (MTJs) are of great interest for non-conventional computing appli-
cations. The Toffoli gate is a universal reversible logic gate, enabling the construction of arbitrary
boolean circuits. Here, we present a proof-of-concept construction of a gadget which encodes the
Toffoli gate's truth table into the ground state of coupled uniaxial nanomagnets that could form the
free layers of perpendicularly magnetized MTJs. This construction has three input bits, three output
bits, and one ancilla bit. We numerically simulate the seven macrospins evolving under the stochastic
Landau-Lifshitz-Gilbert (s-LLG) equation. We investigate the effect of the anisotropy-to-exchange-
coupling strength ratio H<sub>A</sub>/H<sub>ex</sub> on the working of the gadget. We find that for H<sub>A</sub>/H<sub>ex</sub> ≲ 0.93, the
spins evolve to the Toffoli gate truth table configurations under LLG dynamics alone, while higher
H<sub>A</sub>/H<sub>ex</sub> ratios require thermal annealing due to suboptimal metastable states. Under our chosen
annealing procedure, the s-LLG simulation with thermal annealing achieves a 100% success rate up
to H<sub>A</sub>/H<sub>ex</sub> ≈ 3.0. The feasibility of constructing MTJ-free-layer-based Toffoli gates highlights their
potential in designing new types of MTJ-based circuits.

# I. INTRODUCTION

Magnetic tunnel junctions (MTJs) hold significant potential in advancing future computing technologies. Traditionally, MTJs have been employed in the development of long-term storage elements and magnetic random access memory (MRAM). These applications utilize the stable and non-volatile characteristics of MTJs to store data efficiently and reliably [1, 2]. On the other hand, the stochastic behavior and dynamic nature of MTJs have led to novel computing applications beyond data storage [3]. With the rapid advances in machine learning and artificial intelligence, neuromorphic computing has emerged as a promising approach to meet the increasing demand for hardware capable of handling large-scale data processing, resolving the memory bottleneck issue, and having less energy consumption [4, 5]. MTJs have been shown to be a prominent candidate in the development of neuromorphic computing hardware [6–9].

MTJs also provide novel avenues for solving op-
timization problems. Probabilistic Ising machines
built using MTJs have been used to find approxi-
mate solutions to NP-hard problems by encoding the
true solution in the system's ground state [10-12].
Recently, several of the authors of this article con-
ducted numerical studies to benchmark the perfor-
mance of Landau-Lifshitz-Gilbert (LLG) dynamical
systems in solving combinatorial optimization prob-
lems [13]. This work demonstrated that the stochas-
tic and nonlinear dynamics governed by the stochas-

<figure><img src="image_2.png" /><figcaption>FIG. 1. The circuit representation of Toffoli gate (left) and an Ising spins network for the construction of a Toffoli gate (right) [14].</figcaption></figure>

tic LLG equation has the potential to find better so-
lutions to NP problems as system size scales up [13].

Building on these insights, we aim to broaden the
scope of these unconventional computing application
of MTJs. Although some combinatorial optimization
problems, such as 3-SAT, are NP-complete and thus
universal for computation, it can be quite cumber-
some to encode the problem of interest into one of
these Boolean satisfiability problems. It would be
useful to be able to efficiently encode simple tasks,
such as addition, multiplication or other basic oper-
ations into these unconventional dynamical systems.
On the most general level, it suffices to encode a uni-
versal (reversible) logical gate into the ground state
of some coupled MTJs. Here we focus on the Tof-
foli gate, which is a universal reversible gate [15] that
can serve as the building block for arithmetic operations. In the field of magnetism, efforts have been made to develop the Toffoli gate using three interacting classical spins [16], skyrmions [17–18], and spin waves [19]. In this paper, we model the construction of a classical Toffoli gate using seven coupled uniaxial nanomagnets, which could represent the free layers of seven perpendicularly magnetized MTJs. We use the macrospin approximation for each free layer and apply LLG dynamics to simulate their evolution.

In the following sections, we will discuss our methods and results in detail. We begin by describing the mapping of the Toffoli gate’s truth table to the ground state of an Ising Hamiltonian. Next, we outline our approach to simulating the macrospin model using the stochastic Landau-Lifshitz-Gilbert (sLLG) equation and discuss how macrospins can emulate Ising spins. We then detail our method for using coupled macrospins to construct a Toffoli gate. Finally, we present the results of our macrospin Toffoli gate simulations, both at zero temperature and using simulated annealing.

## II. METHODS

### A. Toffoli Gate Construction

The truth table of the Toffoli gate is described by the following Boolean function:

$$ \{a, b, c\} \rightarrow \{a, b, c \oplus (a \land b)\}. \quad (1) $$

That is, the last bit flips conditional on the first two bits being 1. (It can be considered a '2'-controlled-NOT gate.) This gate is universal, as $\{a, b, 1\}$ will result in NAND of the inputs a and b, i.e., $\{a, b, \neg(a \land b)\}$ and the NAND gate is a universal digital logic gate. Since the truth table can be represented as a permutation matrix, the Toffoli gate is also reversible.

This truth table of the Toffoli gate can be mapped to the ground states of a seven-coupled Ising spin system as shown in Chamon et al. [14]. We denote the Ising spins as $s$, which can only take two values $s = \pm 1$. The connection network between Ising spins is depicted in Fig. 1. In this system, we denote the spins that act as inputs $s_a, s_b, s_c$, the spin that acts as an ancilla $s_d$, and the spins that act as the outputs as $s_{a'}, s_{b'}, s_{c'}$. The Hamiltonian for this Toffoli gate construction is given by [14]

$$ H_T = - \sum_{\langle i,j \rangle} J_{ij} s_i s_j - \sum_i h_i s_i, \quad (2) $$

where the $J_{ij}$ is the coupling term between spins $i$ and $j$, and $h_i$ is the local field term applied to each

individual spin. The values for $J_{ij}$ and $h_i$ are shown in Table I.

<table><caption>TABLE I. $J_{ij}$-coupling between the input spins $\{s_a, s_b, s_c\}$, output $\{s_{a'}, s_{b'}, s_{c'}\}$ and the ancilla $s_d$. The diagonal represents the external field $h_i$ applied on the $i^{th}$ spin.</caption><thead><tr><th>s<sub>a</sub></th><th>s<sub>b</sub></th><th>s<sub>c</sub></th><th>s<sub>d</sub></th><th>s<sub>a'</sub></th><th>s<sub>b'</sub></th><th>s<sub>c'</sub></th><th></th></tr></thead><tbody><tr><td>-1</td><td>3</td><td>2</td><td>-4</td><td>1</td><td>0</td><td>-2</td><td>s<sub>a</sub></td></tr><tr><td>3</td><td>3</td><td>-4</td><td>8</td><td>0</td><td>1</td><td>4</td><td>s<sub>b</sub></td></tr><tr><td>2</td><td>-4</td><td>2</td><td>6</td><td>0</td><td>0</td><td>4</td><td>s<sub>c</sub></td></tr><tr><td>-4</td><td>8</td><td>6</td><td>-4</td><td>0</td><td>0</td><td>-6</td><td>s<sub>d</sub></td></tr><tr><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>s<sub>a'</sub></td></tr><tr><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>s<sub>b'</sub></td></tr><tr><td>-2</td><td>4</td><td>4</td><td>-6</td><td>0</td><td>0</td><td>-2</td><td>s<sub>c'</sub></td></tr></tbody></table>

The Ising spins can only take two values: $s = \pm 1$, where $s = +1$ corresponds to logical 1, and $s = -1$ corresponds to logical 0. We denote logical input as $a, b, c$, logical ancilla as $d$, and logical output as $a', b', c'$. For simplicity, the configuration of the system is denoted as a binary string $(a, b, c, d, a', b', c')$, and we label each configuration by the decimal value of this binary string. For example, an input 001, an ancilla 0, and an output 001 form the binary string (0, 0, 1, 0, 0, 0, 1), which corresponds to 17 in decimal representation. With this convention, we plot the energy landscape of the Hamiltonian $H_T$ for different configurations in Fig. 2. The configurations with the lowest energies correspond to the correct input-output relations for the Toffoli gate’s truth table.

Given the Toffoli gate Hamiltonian $H_T$, we model a system of 7 coupled single-domain ferromagnets that represent the 7 Ising spins. The macrospin approximation is applied to each domain such that the magnetization of the spins in the $i^{th}$ domain are represented by one giant magnetic moment $\vec{m}_i = m_i \cdot \hat{m}_i$, with magnitude $m_i$ and unit direction vector $\hat{m}_i$. This giant magnetic moment $\vec{m}_i$ is referred to as a macrospin. We denote the macrospin system as $(\vec{m}_a, \vec{m}_b, \vec{m}_c, \vec{m}_d, \vec{m}_{a'}, \vec{m}_{b'}, \vec{m}_{c'})$. In contrast to Ising spins that are either up or down, a macrospin $\vec{m}$ is defined on the surface of a 3D sphere with two degrees of freedom (in spherical coordinates, $\theta$ and $\phi$), and evolves according to the LLG equation.

To emulate Ising spins using macrospins that follow full classical LLG dynamics, we include a uniaxial magnetic anisotropy along the $\hat{z}$-axis, such that the macrospins favor aligning up or down along the $\hat{z}$-axis. We initialize the system by converting $s = \pm 1$ to $\vec{m} = \pm m \cdot \hat{z}$, aligning the macrospin with the north or south pole based on the sign of the Ising spin. The three macrospins that act as the input
<figure><img src="image_3.png" /><figcaption>FIG. 2. The energy landscape calculated from the Hamil-
tonian $H_T$. The x axis represents the decimal represen-
tation of the binary string (a, b, c, d, a', b', c') and the y
axis is the energy calculated from $H_T$ by converting each
logical bits to its corresponding Ising spins. Each color
represents a set of configurations with the same fixed in-
put values (a, b, c). The blue dots represent the eight
configurations with the lowest energies, which have the
correct input-output relation according to the truth ta-
ble for the Toffoli gate. The unit of energy will be taken
to be $k_B T_{amb}$, where $k_B$ is the Boltzmann constant and
$T_{amb} = 300$ K, in the stochastic macrospin model dis-
cussed in Secs. II C and III B.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. Illustration of macrospin dynamics (left) and its maping to an Ising spin (right).</figcaption></figure>

of the Toffoli gate, ($\vec{m}_a, \vec{m}_b, \vec{m}_c$), are fixed numerically to their initial states throughout the simulation, while the other four macrospins ($\vec{m}_d, \vec{m}_{a'}, \vec{m}_{b'}, \vec{m}_{c'}$) undergo LLG evolution. For readout at the end of the simulation, the sign of the projection of the macrospin's moment on the z-axis (sgn($\hat{m}_z$)) is the binary value of Ising spin, as illustrated in Fig. 3.

B. LLG and s-LLG Equations

To describe the LLG dynamics that each macrospin
follows, we begin by considering the micromagnetic
energy of the system, *E*, which consists of three con-
tributions:

$$
E = E_{\text{ex}} + E_{\text{anis}} + E_{\text{Z}}, \quad (3)
$$

where $E_{\text{ex}}$ is the exchange energy between each pair of moments, $E_{\text{anis}}$ is the anisotropic energy arising from the crystallographic structure of the domain, and $E_{\text{Z}}$ is the Zeeman energy from an external field applied to the moment. The expression for each term is given by:

$$
E_{\text{ex}} = -\mu_0 \sum_{\langle i,j \rangle}^{n} 2J_{ij} H_{\text{ex}} \frac{\vec{m}_i \cdot \vec{m}_j}{m_i + m_j}, \quad (4)
$$

$$
E_{\text{anis}} = -\frac{1}{2}\mu_0 \sum_{i}^{n} H_{Ai} m_i (\hat{e}_{hi} \cdot \hat{m}_i)^2, \quad (5)
$$

$$
E_Z = -\mu_0 \sum_i^n \vec{m}_i \cdot h_i H_{\text{ext}} \hat{e}_{\text{exti}}, \quad (6)
$$

where $\mu_0$ is the vacuum permeability, $J_{ij}$ are dimensionless coupling terms characterizing the exchange interaction between magnetic moments, and $H_{\text{ex}}$ (with units of field) characterizes the general scale of the exchange coupling strength. The summation over $\langle i, j \rangle$ is over pairs of connected spins. $H_{Ai}$ characterizes the magnetic anisotropy of moment $i$, and $\hat{e}_{hi}$ is the unit vector representing the preferred crystallographic direction of moment $i$, which we set $\hat{e}_{hi} = \hat{z}$ to facilitate readout along this axis. $h_i$ are dimensionless terms that characterize the external magnetic field applied to each magnetic moment, and $H_{\text{ext}}$ (with units of field) characterize the general scale of the external field strength. $\hat{e}_{\text{exti}}$ is the unit vector in the direction of the applied external field. We set $\hat{e}_{\text{exti}} = \hat{z}$, so that this field is along the same axis as the magnetic anisotropy.

The $i^{th}$ magnetic moment experiences an effective field $\vec{H}_{\text{eff},i}$, which can be obtained by computing the derivative of the micromagnetic energy $E$ with respect to the magnetic moment $\vec{m}_i$:

$$
\vec{H}_{\text{eff},i} = - \frac{1}{\mu_0} \frac{\partial E}{\partial \vec{m}_i}. \qquad (7)
$$

The expression for $\vec{H}_{\text{eff},i}$ is substituted into LLG
equation, which determines the dynamical evolution
of each magnetic moment at each time $t$:

$$
\frac{d\vec{m}_i}{dt} = \gamma' \vec{m}_i \times \vec{H}_{\text{eff},i} - \frac{\alpha_i}{m_i} \vec{m}_i \times \frac{d\vec{m}_i}{dt}, \quad (8)
$$
where $\gamma' = \mu_0\gamma$, $\gamma$ is the gyromagnetic ratio. $\alpha_i$ is the damping constant of the $i^{th}$ moment.

To account for the stochastic nature of the system under thermal fluctuations, a Langevin field term $\vec{H}_{th,i}$ is added to the effective field $\vec{H}_{eff,i}$:

$$ \vec{H}_{eff,i} \rightarrow \vec{H}_{eff,i} + \vec{H}_{th,i}. \quad (9) $$

The Langevin fields at times $t$ and $t'$ follow the statistical properties described in Refs. [20–22]:

$$ \langle \vec{H}_{th,i} \rangle = 0, \quad (10) $$

$$ \langle \vec{H}_{th,i,m}(t) \vec{H}_{th,i,n}(t') \rangle = C_i \delta_{mn} \delta(t-t'), \quad (11) $$

$$ C_i = \frac{2\alpha_i k_B T}{\mu_0^2 m_i \gamma (1 + \alpha_i^2)}, \quad (12) $$

where $m$ and $n$ represent the $x, y, z$ components of the random field, and $\delta_{mn}$ and $\delta(t-t')$ are the Kronecker and Dirac delta functions respectively. $t$ and $t'$ denote different time instances, $k_B$ is the Boltzmann constant, and $T$ is the temperature. By integrating Eq. (11) over the discretized time step $\Delta t$, we can derive the standard deviation of $\vec{H}_{th,i}$ in the $x, y, z$ directions:

$$ \sigma_i = \frac{1}{\mu_0} \sqrt{\frac{2\alpha_i k_B T}{m_i \gamma (1 + \alpha_i^2) \Delta t}}. \quad (13) $$

For each component $x, y, z$ of $\vec{H}_{th,i}$, we sample a value from a Gaussian distribution with a mean $\mu = 0$ and standard deviation given by $\sigma_i$ from Eq. (13).

For the numerical integration, we use the Stratonovich interpretation of the stochastic process and implement Heun's method as our integration scheme. In the study by Ament et al. [23], Heun's integration scheme converges to the Stratonovich solution and preserves the norm of $\hat{m}_i$ for small time steps when simulating s-LLG dynamics.

### C. Implementation of Toffoli Gate Hamiltonian

In our simulation, we assume identical parameters for all $n$ macrospins, allowing us to omit the index $i$ for the following quantities: $H_{Ai} = H_A$, $\alpha_i = \alpha$, and $m_i = m$. This approach models an array of perfectly identical macrospins. The simulation can then be simplified to focus solely on the evolution of the unit magnetic moment $\hat{m}_i$, as shown by the blue $\hat{m}$ trajectories in Fig. 3. The energy terms can thus be

simplified as follows:

$$ E_{\text{ex}} = -\mu_0 m H_{\text{ex}} \sum_{\langle i,j \rangle}^{n} J_{ij} (\hat{m}_i \cdot \hat{m}_j), \quad (14) $$

$$ E_{\text{anis}} = -\frac{1}{2} \mu_0 m H_A \sum_{i}^{n} (\hat{m}_i \cdot \hat{z})^2, \quad (15) $$

$$ E_Z = -\mu_0 m H_{\text{ext}} \sum_{i}^{n} h_i (\hat{m}_i \cdot \hat{z}). \quad (16) $$

The Toffoli gate Hamiltonian $H_T$ can be directly mapped to the system's micromagnetic energy terms $E_{ex}$ and $E_Z$. Specifically, the $J_{ij}$ terms in $H_T$ in Eq. (2) correspond to the $J_{ij}$ terms in $E_{ex}$ in Eq. (14), and the $h_i$ terms in $H_T$ correspond to the $h_i$ terms in $E_Z$ in Eq. (16). To maintain the relative strength between the coupling terms and the local field terms as given by $H_T$, we set the general scale of the exchange coupling strength $H_{ex}$ equal to the local field strength $H_{ext}$, i.e., $H_{ex} = H_{ext}$ for all trials.

The ratio of anisotropy $H_A$ to coupling strength $H_{ex}$ influences the quality of the final output. High anisotropy causes spins to align more strictly along the north and south poles, emulating binary Ising spins. However, high anisotropy also makes it more difficult for spins to overcome the energy barrier to flip, increasing the likelihood that the system will be trapped in configurations that are not consistent with the truth table of the Toffoli gate. Conversely, if $H_A$ is relatively small, the spins can flip more easily, but this can result in the spins not aligning well with the $z$-axis, complicating the read-out process. In addition, thermal effects play a role in reaching the correct ground state when simulated annealing is incorporated. High temperatures produce thermal fluctuations that allow spins to overcome energy barriers and eventually find the ground state through annealing.

Considering the many possible contributing factors, we simulate the zero-temperature dynamics and the dynamics with thermal annealing separately. First, we investigate the ability of pure (i.e., zero temperature) LLG dynamics to find the correct configurations given by the truth table as a function of $H_A/H_{ex}$, with $H_{ex} = H_{ext}$ fixed. The values of the parameters for our simulation are shown in Table II. The values are realistic for the free layers of perpendicularly magnetized magnetic tunnel junction nanopillars [1].

In addition, we study the ability of combined stochastic LLG dynamics with simulated annealing to reach the true ground state as a function of $H_A/H_{ex}$. We have empirically chosen an annealing schedule that decreases from $T = 300$ to 0K over 150 evenly
spaced steps. We set the exchange field $H_{ex}$ such that the exchange energy $E_{ex}$ is the order $k_B T_{amb}$, with $T_{amb} = 300$ K. This ensures the system has sufficient thermal energy to overcome energy barriers as it evolves. At each temperature, we let the system evolve with $10^5$ iterations, which corresponds to $1\mu s$ in physical units. This time is longer than a typical spin relaxation time. For instance, considering spin relaxation in the exchange field, $t_{relax} \approx 2\pi/(\alpha\mu_0\gamma H_{ex}) \approx 0.1\mu s$, which is one order magnitude smaller than $1\mu s$. Therefore, our system will be in or very close to thermal equilibrium with the environment at the end of each temperature step.

<table><caption>TABLE II. Values of the parameters chosen for simulation.</caption><tr><td>Damping Constant</td><td>α</td><td>0.01</td></tr><tr><td>Magnetic Moment</td><td>m</td><td>8 × 10<sup>-19</sup> J/T</td></tr><tr><td>Ambient Temperature</td><td>T<sub>amb</sub></td><td>300K</td></tr><tr><td>Exchange Field</td><td>H<sub>ex</sub></td><td>4.12 × 10<sup>3</sup> A/m</td></tr><tr><td>Simulation Time Step</td><td>Δt</td><td>0.01 ns</td></tr><tr><td>Iterations at Each Temperature</td><td></td><td>1 × 10<sup>5</sup></td></tr><tr><td>Iterations at 0K (LLG)</td><td></td><td>5 × 10<sup>5</sup></td></tr></table>

### III. RESULTS

#### A. Zero temperature simulations

We first simulate the construction of a Toffoli Gate using LLG dynamics at zero temperature ($T = 0$ K). For a given initial configuration $(a, b, c, d, a', b', c')$, we begin by initializing the macrospins to the corresponding values $(\hat{m}_a, \hat{m}_b, \hat{m}_c, \hat{m}_d, \hat{m}_{a'}, \hat{m}_{b'}, \hat{m}_{c'})$. The input bits $(a, b, c)$ are converted into macrospins $(\hat{m}_a, \hat{m}_b, \hat{m}_c)$ by aligning them strictly along the $\hat{z}$ axis, with bit 0 corresponding to $-\hat{z}$ and bit 1 to $+\hat{z}$. Then these three macrospins are numerically fixed throughout the simulation. We initialize the output bits $(d, a', b', c')$ as $(\hat{m}_d, \hat{m}_{a'}, \hat{m}_{b'}, \hat{m}_{c'})$ by first aligning these vectors along the $z$ axis. Then we add a small angular deviation of 0.01 rad from $\hat{z}$ to prevent the spins from experiencing zero initial torque (i.e., to avoid the cross-product terms in Eq. (8) being zero).

For each Toffoli Gate input $(a, b, c)$, there are 16 different possible initial configurations of $(d, a', b', c')$. For a fixed input $(a, b, c)$, we run 16 trials with different initial conditions of $(d, a', b', c')$, allowing the system to evolve to its final state. We record the final state and count the number of trials that result in an output $(a', b', c')$ consistent with the Toffoli gate truth table and use this to compute the success rate.

<figure><img src="image_5.png" /><figcaption>FIG. 4. Success rate of the Toffoli gate reaching the correct output as a function of $H_A/H_{ex}$ ratio, under deterministic (zero temperature) LLG dynamics. The results for each input $(a, b, c)$ are color-coded in the plot. The graph appears to have fewer colors than listed due to multiple overlapping lines stacked on top of each other. The black dashed line at $H_A/H_{ex} \approx 0.93$ represents the largest $H_A/H_{ex}$ in our simulation that still achieves a 100% success rate for all trials.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Average $\hat{z}$ projection of output spins $\langle|\hat{m}_z|\rangle$ as a function of $H_A/H_{ex}$ ratio, under deterministic (zero temperature) LLG dynamics. The results for each input configuration $(a, b, c)$ are color-coded in the plot. The graph appears to have fewer colors than listed due to multiple overlapping lines stacked on top of each other. The black dashed line at $H_A/H_{ex} \approx 0.93$ represents the largest $H_A/H_{ex}$ in our simulation that still achieves a 100% success rate for all trials.</figcaption></figure>

The process is repeated for different $H_A/H_{ex}$ ratios, and the result is shown in Fig. 4. We observe a 100%
success rate of reaching the correct configuration for $H_A/H_{ex} \lesssim 0.93$.

To assess how Ising-like the spins are at the end of the simulation, we record the absolute magnitude of the final $\hat{z}$ projection $|\hat{m}_z| = |\hat{m} \cdot \hat{z}|$ for each output configuration $(\hat{m}_d, \hat{m}_{a'}, \hat{m}_{b'}, \hat{m}_{c'})$. There are 16 different trials for each input $(a, b, c)$, we take the average of $|\hat{m}_z|$ across both the four output spins and the 16 trials to obtain a single value. The results are shown in Fig. 5. Since we set the external field applied locally to each macrospin as $\hat{e}_{ext} = \hat{z}$, an Ising-like output is still observed even when $H_A = 0$, with the smallest $\langle|\hat{m}_z|\rangle = 0.72$. We also observe that as $H_A/H_{ex}$ increases, $\langle|\hat{m}_z|\rangle$ approaches 1, which indicates that the spins are more Ising like as the anisotropy increases.

## B. Simulation with thermal annealing (s-LLG)

In this section, we simulate the Toffoli gate construction with simulated annealing. We fix the exchange coupling strength such that $E_{ex}$ between spins is on the order of $k_B T_{amb}$, where $T_{amb}$ is the ambient temperature of 300K. For each Toffoli Gate input $(a, b, c)$, we start with the configuration $(a, b, c, 0, 0, 0, 0)$ and perform an annealing process from 300 to 0K over 150 evenly spaced steps. Since we start at high temperature, the system is able to visit different magnetic configurations. Therefore, it is not as necessary to initialize the system with different configurations $(d, a', b', c')$, as in the 0K case. After the annealing process, we record the final state and repeat the simulations for 500 independent trials. We do the same simulations for different ratios $H_A/H_{ex}$. The results are shown in Figs. 6 and 7.

With the assistance of simulated annealing, we observe that within the statistics of 500 trials, the system is able to achieve a 100% probability of finding the correct output for $H_A/H_{ex}$ up to $H_A/H_{ex} \approx 3.0$, which is higher than the limit of $H_A/H_{ex} \approx 0.93$ observed with zero-temperature pure LLG dynamics. Further, at $H_A/H_{ex} \approx 3.0$, the final configurations are all Ising-like, with all spins aligned along the $z$ axis as shown in Fig. 7.

## C. Conclusion

We have demonstrated a proof-of-concept Toffoli gate construction using macrospins with a uniaxial magnetic anisotropy that could be formed by the free magnetic layers of MTJs. Specifically, we investigated the influence of the ratio of the uniaxial magnetic anisotropy to the exchange coupling, $H_A/H_{ex}$,

<figure><img src="image_7.png" /><figcaption>FIG. 6. Probability of reaching the correct configurations at different $H_A/H_{ex}$ ratios for Toffoli gate construction with s-LLG dynamics and simulated annealing. The results for each $(a, b, c)$ are color-coded in the plot. For each $(a, b, c)$, the simulation is repeated for 500 times. Within the statistics of 500 trials, values above $H_A/H_{ex} \approx 3.0$ fails to reach a 100% probability of producing the correct output, as shown by the blue dashed line.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Average $\hat{z}$ projection of output spins, $\langle|\hat{m}_z|\rangle$, as a function of $H_A/H_{ex}$ ratios for Toffoli gate construction with s-LLG dynamics and simulated annealing. The average is taken across 500 trials. The results for each $(a, b, c)$ are color-coded in the plot. The blue dashed line at $H_A/H_{ex} \approx 3.0$ represents the largest $H_A/H_{ex}$ allowed to achieve a 100% probability of reaching the correct output with simulated annealing.</figcaption></figure>

on the spin dynamics. Remarkably, the zero temperature LLG dynamics always ends in configurations consistent with the Toffoli gate's truth table provided that $H_A/H_{ex} \lesssim 0.93$. Introducing thermal annealing
ing allows to push this ratio to higher values, up to $H_A/H_{ex} \simeq 3.0$, while keeping a 100% success rate.

Our modeling thus shows that a universal reversible gate can be constructed from interacting classical macrospins, highlighting a promising avenue for research with arrays of interacting MTJs [24–26]. Although our simulation has shown the feasibility of this approach, further research is needed. Specifically, this work explored the feasibility of constructing a single Toffoli gate. Future theoretical work could investigate the feasibility of building circuits using multiple Toffoli gates for real-world computational tasks. These would address the questions of methods, properties, and fidelity of connected Toffoli gates and, for example, how the time to solution

scales with the problem size.

Acknowledgement

The research at NYU was supported by the Office of Naval Research (ONR) under Award No. N00014-23-1-2771 and in part by the National Science Foundation (NSF) under Award DMR-2105114. Research at the University of Lorraine was supported by “MAT-PULSE,” part of the French PIA project “Lorraine Université d’Excellence” reference ANR-15-IDEX-04-LUE. This work was also supported in part through the NYU IT High Performance Computing resources, services, and staff expertise.

[1] A. D. Kent and D. C. Worledge, A new spin on magnetic memories, Nature Nanotechnology 10, 187 (2015).

[2] J. Ender, S. Fiorentini, R. L. De Orio, W. Goes, V. Sverdlov, and S. Selberherr, Emerging CMOS Compatible Magnetic Memories and Logic, IEEE Journal of the Electron Devices Society 9, 456 (2021).

[3] G. Finocchio, J. A. C. Incorvia, J. S. Friedman, Q. Yang, A. Giordano, J. Grollier, H. Yang, F. Ciubotaru, A. V. Chumak, A. J. Naeemi, S. D. Cotofana, R. Tomasello, C. Panagopoulos, M. Carpentieri, P. Lin, G. Pan, J. J. Yang, A. Todri-Sanial, G. Boschetto, K. Makasheva, V. K. Sangwan, A. R. Trivedi, M. C. Hersam, K. Y. Camsari, P. L. McMahon, S. Datta, B. Koiller, G. H. Aguilar, G. P. Temporão, D. R. Rodrigues, S. Sunada, K. Everschor-Sitte, K. Tatsumura, H. Goto, V. Pulitafito, J. Åkerman, H. Takesue, M. D. Ventra, Y. V. Pershin, S. Mukhopadhyay, K. Roy, I.-T. Wang, W. Kang, Y. Zhu, B. K. Kaushik, J. Hasler, S. Ganguly, A. W. Ghosh, W. Levy, V. Roychowdhury, and S. Bandyopadhyay, Roadmap for unconventional computing with nanotechnology, Nano Futures 8, 012001 (2024).

[4] C. D. Schuman, S. R. Kulkarni, M. Parsa, J. P. Mitchell, P. Date, and B. Kay, Opportunities for neuromorphic computing algorithms and applications, Nature Computational Science 2, 10 (2022).

[5] W. Zhang, B. Gao, J. Tang, P. Yao, S. Yu, M.-F. Chang, H.-J. Yoo, H. Qian, and H. Wu, Neuroinspired computing chips, Nature Electronics 3, 371 (2020).

[6] J. Torrejon, M. Riou, F. A. Araujo, S. Tsunegi, G. Khalsa, D. Querlioz, P. Bortolotti, V. Cros, K. Yakushiji, A. Fukushima, H. Kubota, S. Yuasa, M. D. Stiles, and J. Grollier, Neuromorphic computing with nanoscale spintronic oscillators, Nature 547, 428 (2017).

[7] M. Romera, P. Talatchian, S. Tsunegi, F. Abreu Araujo, V. Cros, P. Bortolotti, J. Trastoy, K. Yakushiji, A. Fukushima, H. Kubota, S. Yuasa, M. Ernoult, D. Vodenicarevic, T. Hirtzlin, N. Lo-catelli, D. Querlioz, and J. Grollier, Vowel recognition with four coupled spin-torque nano-oscillators, Nature 563, 230 (2018).

[8] J. Grollier, D. Querlioz, K. Y. Camsari, K. Everschor-Sitte, S. Fukami, and M. D. Stiles, Neuromorphic spintronics, Nature Electronics 3, 360 (2020).

[9] D. Marković, A. Mizrahi, D. Querlioz, and J. Grollier, Physics for neuromorphic computing, Nature Reviews Physics 2, 499 (2020).

[10] K. Y. Camsari, R. Faria, B. M. Sutton, and S. Datta, Stochastic p-Bits for Invertible Logic, Physical Review X 7, 031014 (2017).

[11] W. A. Borders, A. Z. Pervaiz, S. Fukami, K. Y. Camsari, H. Ohno, and S. Datta, Integer factorization using stochastic magnetic tunnel junctions, Nature 573, 390 (2019).

[12] J. Si, S. Yang, Y. Cen, J. Chen, Y. Huang, Z. Yao, D.-J. Kim, K. Cai, J. Yoo, X. Fong, and H. Yang, Energy-efficient superparamagnetic Ising machine and its application to traveling salesman problems, Nature Communications 15, 3457 (2024).

[13] D. Chen, A. D. Kent, D. Sels, and F. Morone, Solving combinatorial optimization problems through stochastic Landau-Lifshitz-Gilbert dynamical systems (2024), arXiv:2407.00530.

[14] C. Chamon, E. R. Mucciolo, A. E. Ruckenstein, and Z.-C. Yang, Quantum vertex model for reversible classical computing, Nature Communications 8, 15303 (2017).

[15] T. Toffoli, Reversible computing, in *Automata, Languages and Programming*, edited by J. de Bakker and J. van Leeuwen (Springer, Berlin, Heidelberg, 1980) pp. 632–644.
[16] D. Nuzzi, L. Banchi, R. Vaia, E. Compagno, A. Cucoli, P. Verrucchi, and S. Bose, Full-magnetic implementation of a classical toffoli gate, Phys. Rev. Res. 6, 013308 (2024).

[17] J. I. Costilla, J. W. Alegre, A. Talledo, and B. R. Pujada, Implementation of the Toffoli and Peres reversible logic gates using magnetic skyrmions in operational gates, Journal of Applied Physics 134, 013903 (2023).

[18] M. Chauwin, X. Hu, F. Garcia-Sanchez, N. Betrabet, A. Paler, C. Moutafis, and J. S. Friedman, Skyrmion Logic System for Large-Scale Reversible Computation, Physical Review Applied 12, 064053 (2019).

[19] M. Balynskiy, H. Chiang, D. Gutierrez, A. Kozhevnikov, Y. Filimonov, and A. Khitun, Reversible magnetic logic gates based on spin wave interference, Journal of Applied Physics 123, 144501 (2018).

[20] W. F. Brown, Thermal Fluctuations of a Single-Domain Particle, Physical Review 130, 1677 (1963).

[21] D. V. Berkov, Magnetization Dynamics Including Thermal Fluctuations: Basic Phenomenology, Fast Remagnetization Processes and Transitions Over High-energy Barriers, in Handbook of Magnetism and Advanced Magnetic Materials (John Wiley & Sons, Ltd, 2007).

[22] J. L. García-Palacios and F. J. Lázaro, Langevin-dynamics study of the dynamical properties of small magnetic particles, Phys. Rev. B 58, 14937 (1998).

[23] S. Ament, N. Rangarajan, A. Parthasarathy, and S. Rakheja, Solving the stochastic landau-lifshitz-gilbert-slonczewski equation for monodomain nano-magnets: A survey and analysis of numerical techniques (2017), arXiv:1607.04596

[24] P. Talatchian, M. W. Daniels, A. Madhavan, M. R. Pufall, E. Jué, W. H. Rippard, J. J. McClelland, and M. D. Stiles, Mutual control of stochastic switching for two electrically coupled superparamagnetic tunnel junctions, Physical Review B 104, 054427 (2021).

[25] N.-T. Phan, L. Soumah, A. Sidi El Valli, L. Hutin, L. Anghel, U. Ebels, and P. Talatchian, Electrical Coupling of Perpendicular Superparamagnetic Tunnel Junctions for Probabilistic Computing, in Proceedings of the 17th ACM International Symposium on Nanoscale Architectures (ACM, Virtual OR USA, 2022) pp. 1-6.

[26] L. Schnitzspan, M. Kläui, and G. Jakob, Electrical coupling of superparamagnetic tunnel junctions mediated by spin-transfer-torques, Applied Physics Letters 123, 232403 (2023).