Logical gates embedding in Artificial Spin Ice

The realization and study of arrays of interacting magnetic nanoislands, such as artificial spin ices, have reached mature levels of control that allow design and demonstration of exotic, collective behaviors not seen in natural materials. Advances in the direct manipulation of their local, binary moments also suggest a use as nanopatterned, interacting memory media, for computation {\it within} a magnetic memory. Recent experimental work has demonstrated the possibility of building logic gates from clusters of interacting magnetic domains, and yet the possibility of large scale integration of such gates can prove problematic even at the theoretical level. Here we introduce theoretically complete sets of logical gates, in principle realizable in an experiment, and we study the feasibility of their integration into tree-like circuits. By evaluating the fidelity control parameter between their collective behavior and their expected logic functionality we determine conditions for integration. Also, we test our numerical results against the presence of disorder in the couplings, showing that the design gate structure is robust to small coupling perturbations, and thus possibly to small imperfections in the fabrication of the islands.


Introduction
The last fourteen years have seen the use of interacting [1,2] magnetic nanostructures [3] patterned in different geometries in so called artificial spin ices, to realize a wealth of different emergent behaviors often not found in natural magnets. The level of control afforded by these materials has allowed the study of frustration and residual entropy, within a broad range of exotic phenomena and potential applications [4,5,6,4,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]. Recently, Gartside et al. and Wang et al. have demonstrated the fine, local manipulation of islands magnetization [23,24,25] to reliably write the bites of artificial spin ice [1,20] materials. Moreover, still unpublished results [28,29] demonstrate local activation of the islands kinetics via photo-induced heating. The use of artificial spin ices as collective rather than individual forms of memory has also been explored [30,31].
This confluence of results opens now intriguing perspectives toward new platforms for computation within a memory medium. In current magnetic memory storage, it is important that each magnetic domain is not influenced by the other domains to prevent data corruption. However, an interacting memory could allow for computation performed within the memory [32,33,34]. In a network of interacting memory bits there would be no obvious directionality of flow from input to output within the logic gates composing it. It has been proposed that, in the absence of a logic unit, information overhead (the capability of compressing information in the collective state of the network of magnets), intrinsic parallelism, and functional polymorphism would be natural [35].
In a network of interacting, binary spins, some could be kept fixed in time and considered inputs via external fields or via nearby larger islands used as biases. Other islands could be considered outputs or results intermediate computations, to be read after the system has relaxed in a low-energy state; this is some other fixed point, which would correspond to the result of a computation. Indeed, Boolean gates of interacting, magnetic nanoislands have been already pioneered [36,37,38] in the context of nonstandard logic and more recently within artificial spin ice materials [39,40]. However, the crucial issue neglected in most experimental work is the reliability of their integration in realistic circuits, which can be problematic [41,42,39].
In this work we explore theoretically the feasibility of integrating nanomagnetic gates in functional circuits by simulating thermal annealing of tree-like circuits made of magnetic nanoislands. Previous work in this respect has been done via the implementation of boolean gates directly in the couplings [26]. However, in this work we consider the possibility of reprogramming the gates a posteriori using intermediate islands exerting a biasing field (see below). This prompts us to perform a new analysis in order to go beyond previous works on the Husimi cactus [27]. We will also consider localized, rather than homogeneous, annealing that proceeds along the circuit, now a possibility via photo-activation [28]. The purpose of this paper is to verify the fidelity control parameter between the convergence of such tree-like interacting systems and the logical functionality of the corresponding logic circuits.

Gates
It is not difficult to devise logical gates from interacting nanoislands whose magnetization can be modeled as a classical, binary Ising spin. For definiteness, consider the NAND and NOR gates (which are universal gates) realized via triangular plaquettes of antiferromagnetically interacting classical Ising spins in Fig. 1 which can be (and have been, though in a different context [43]) fabricated at the nanoscale with magnetic moments perpendicular to the array. We imagine that red spins can be set as inputs and held fixed during computation, while the green spin is the output to be read. The blue spins provide a bias, that is a local field h on the output that eliminates any indeterminacy of the output that results from frustration. Its orientation determines whether the gate is (N)AND or (N)OR, and allows for reconfigurable circuitry. In practice these gates can be realized through a nanoisland of higher coercivity, such that its moment does not change during computation. Switching the magnetization of such an island leads to changing between (N)AND and (N)OR, thus reprogramming the Schematics of a circuit associated with a local magnetization on the triangles (dashed lines represent the interactions, double-arrows the fluctuating spins, blue/red dots the orientation of bias). We highlighted what we called horizontal interactions. Figure 2. A magnetic NAND gate (left) and NOR gate (right) built from antiferromagnetically interacting moments, and their truth table in Boolean algebra and their realization through antiferromagnetic spins. In red we have the input spins, in green the output. The blue spin is needed to bias frustration and changes the plaquette from a NAND to a NOR, making the gate re-programmable.

circuit.
One can then integrate gates into circuits. As an example, Fig. 2 shows integration to produce a 3-bits circuit meant to select Fibonacci numbers.
The goal of this paper is to implement a logic gate structure, such as the one of Fig. 1, in spin ice systems. In Fig. 2, we have shown as an example the case of antiferromagnetic bits; in reality, much more freedom is allowed in choosing the couplings among spins, and therefore in engineering gates, as we show below. Spins can be perpendicular to the array or in plane, and in both cases, have been realized in practice [43]. When spins are in plane, because of the anisotropy of the dipolar interaction, the different relative orientation of the magnetic island can be translated in different coupling constant, ferro-or antiferromagnetic. Indeed it is even possible to remove the coupling between the inputs of a gate setting the moments perpendicularly. As an example, in Fig.3 we represent a Fibonacci series solver using such mapping. In fact, any (N)AND, (N)OR gates can be realized in practice (see Supp. Mat. for a more general analysis).
The gate functionality for three spins can be obtained by studying a simple spin Hamiltonian of the form in which we assume that the spins σ 1 and σ 2 are inputs, and σ 3 is our output and free to fluctuate at a certain temperature T . However, it is not hard to see that if 2|J| > h the ground state of the system are those of a logic gates (N)AND and (N)OR, where the negation N depends on the ferro-or antiferromagnetic coupling, and the type of gates on the biasing field h. In Fig. 4 (f), we show the gate functionality as a function of the parameters J and h, keeping in mind that we assume the underlying structure to be based on magnetic nanoislands and artificial spin ices.

Integration
While gates are easy to conceive in a system of interacting binary variables, there are, however, various problematic aspects in going from a single gate to a circuits. One is the tolerance to errors in the input states, which we discuss in Supp. Mat. A more relevant problem, however, is degeneracy induced by frustration and long-range interactions and/or disorder. Note that gates were introduced under the assumption that the input is held fixed while the output is computed. That cannot be done for gates inside a circuit, where one gate's input is another gate's output. Also, in the absence of frustration or long-range interactions, if the strength of couplings is the same along the hierarchy, then some spins at intermediate stage can arrange itself without changing the output (see Sup. Mat.). While this might seem irrelevant, logical fidelity not only of output, but of each constitutive gate is necessary to go beyond a Turing approach toward neuromorphic computing in a network circuit. This problem can be solved by modulating the couplings J among islands and make them scale with the layer of the tree, which we label by k. One possibility is to choose J k+1 = J k /(2 + ) for some > 0 and h k+1 = h k /2 where h is the local field from the biasing island which selects gate functionality. For the model without horizontal (input) spin interactions, we choose such that |J k /(2 + )| < |h k | and We perform numerical simulations of randomly chosen gates via a Glauber [46] spin dynamics with exponential annealing. The probability of a spin flip leading to an energy change ∆E is given by p(k) = (1 + e ∆E/T k (t) ) −1 , where T k is the temperature at the layer k. We perform both uniform annealing, corresponding to T k (t) = T 0 λ t , with λ < 1 a measure of how slow the annealing is (slower for λ closer to 1), but also a layer dependent temperature to simulate light heating [28,29] orT k (t) = T 0 λ t /2 k : the lower layer (closer to inputs) will cool first, and the top layer last.
As a control parameter for the fidelity of the circuit we introduce a control parameter to measure the difference between the spin system and the equivalent logic circuit, not only for final output, but for all the intermediate layers. For each gate, we consider the vector of spin orientations L corresponding to expected logical functionality, and the one obtained from the interacting spin system, S. Using the vector L, we can define the gate fidelity control parameter as the quantity where N is the total number of gates. Full functionality corresponds to Q = 1, and a system in which all gates fail has Q = 0. Completely random (high temperature) systems have Q = 0.5, as an uncorrelated gate has probability 0.5 of being in the right output state. We call this overlap total. Later we will also use an output overlap, where only the outputs are considered and not the intermediate spins.
We consider two kinds of trees. One kind lacks horizontal interactions (i.e. those between the inputs of gates) annealed homogeneously (HT) and non-homogeneously (NT). The other has horizontal interaction, and also is annealed homogeneously (HTh) and non-homogeneously (NTh). Trees have n = 12 layers of gates, corresponding to 2048 total gates, 2048 floating spins, and input of 2048 bits (e.g. 4096 total spins between inputs and thermally activated spins). On each of the four case we perform 100 simulations corresponding to 100 different circuits (that is, assigning randomly reprogramming biases in each gate) with 100 different random inputs. Fig. 5 shows the average gate overlap for the four cases. In absence of horizontal interactions, the system converges to the proper logical output with total overlap one. This demonstrates gates integration for deterministic computing.
Instead, in the case of horizontal interactions, Q does not converge to one, and the system performs at only about 85% accuracy. Furthermore, as Fig. 6 (top) shows, a noticeable difference can be observed between the curve of homogeneous and nonhomogeneous annealing. The non-homogeneous case converges considerably faster. Also, Fig. 6 (bottom) shows that a non-homogeneously annealed tree with horizontal interaction converges to slightly higher values of overlap for slower annealing, unlike the homogeneous annealing.
It is important to note that the measure Q is a quite restrictive one. One might argue that the result of the computation might still be correct when the intermediate computations contains faulty gates. One can in fact introduce another measure of overlap Q(output) = s out l out + (1 − s out )(1 − l out ), where l out is the expected final result of the computation and s out is the fluctuating (output) spin in our Monte Carlo Figure 5. Behavior of the fidelity control parameter as a function of time when the system is annealed. We denote with mark x curves belonging to exact trees and o curves with horizontal interactions.The curves are obtained after averaging over 100 samples of random gates with n = 10 layers (2048 thermally activated spins and 2048 inputs) random gates. We observe that while for a tree the curves relax to Q = 1, and the trees with horizontal interactions relax to an approximate value of Q ≈ 0.85, thus presenting a macroscopic number of defects. Here 1 time unit is N Monte Carlo flip attempts, with N the number of spins. We also observe that in the case with homogeneous temperature, the fidelity control parameter relaxes a bit slower compared to the non-homogeneous case. In the legend, λ is the annealing rate, e.g. T k = λT k−1 , and thus for λ = 1 no annealing occurs.
Legend: homogeneous temperature on a tree (HT), non-homogeneous temperature on a tree (NT), non-trees with homogeneous (HTh) and non-homogeneous temperature (NTh). In the legend, N stands for nonhomogeneous while H for homogeneous, and T for temperature; h stands for horizontal interactions added. A video of the annealing is provided in [47].
simulations. It is interesting to note that we observe numerically Q(output) > Q . Figure 6, bottom, shows that while loop-containing tree circuits overlap with predicted  As we see, in both cases we obtain that the ovelap is reduced, and already for n = 9 and σ = 1 the fidelity control parameter reaches the completely random value Q = 0.5.
logic functionality at only ∼ 85%, the overlap of the output is still very close to 1, meaning that that internal defects seem to have a finite correlation length (e.g. do not propagate to the output). While loops are detrimental for the general efficiency of the embedding, they do not seem to affect the output. Finally, we tested robustness of computation against quenched disorder. This can come from fabrication, but also from the next nearest neighbor interactions. The strength and signs of these interactions depend on the relative orientation of the islands, and thus not only on the topology but the overall orientation of the circuit. In order to test the role of spurious interactions, we can introduce noise in the couplings and see how this effects the degree of computation. For this purpose, to the overall coupling matrix we add a Gaussian independent and identically distributed noise J σ ij = N (0, σ), e.g. random noise with zero mean and variance σ. In Fig. 7 we show the fidelity control parameter as a function of the noise strength σ, averaged over different realization of the noise. The results thus show that noise can dramatically affect computation in these logical circuits.
We conclude that integration of gates of floating spins into tree-like circuits without loops due to horizontal interactions leads to the theoretical possibility of deterministic computations, whereas loops imply circa 15% faultiness. There is, however, fault tolerance in the output. We interpret this result as an effect of the non-triviality of the energy landscape when the loops are present, and thus to local minima in which the system gets trapped.

Faults, Defects, and Kibble-Zurek
One of the motivations for computing within the memory lies in overcoming the Turing paradigm. That would entail in general loop-containing networks (which might be used to enforce constraints on computation), and also functional polymorphism [33], which can blur the difference between inputs and outputs. It is therefore relevant to investigate the origin of the internal faultiness seen in interaction graphs containing loops. A possible explanation lies in the Kibble-Zurek (KZ) mechanism [49,50], which describes the formation of a nonzero density of defects (excitations about the ground state) when a system crosses a critical temperature during an annealing. The second possibility is a glass transition in the spirit of the random field Ising model, in which the local external fields (our biases) are random. This scenario has been studied in [51] for the case of Cayley trees. There, it was noticed that the random field Ising model has a second order phase transition. The presence of a phase transition is therefore important and worth exploring. As we anneal the system and we cross a transition temperature, defects form at transition either via a coarsening or a Kibble-Zurek mechanism. These defects (and the rate at which they can be reabsorbed) prevent reaching the designed ground state, or in our language, a correct result of the computation.
In order to identify a possible second order phase transition we study the Yang-Lee zeros of the partition function [52], which accumulate toward the real axes of the complex plane (in the thermodynamic limit) when a transition is present. Because we are dealing with trees, it is possible to obtain a recursive formula for the partition function by exploiting an old trick, the star triangle relationship, shown in Fig. 8 (See Appendix). Following for instance Baxter [48], we see that each local triangle can be turned into a local tree via the introduction of a virtual spin in the gates, with coupling constants L 12 , L 3 given by and v 1 = v 2 = tanh(J 3 ), v 3 = tanh(J 12 ). Using these parameters, we obtain a recursion relation for the partition function from which zeros can be computed as a function of the length n of the tree (herẽ h = h/T,L = L/T ). We can thus understand the difference in the annealings of Fig. 5 as follows. Whenever J 12 = 0 , Fig. 9 shows that the zeros in the complex fugacity plane y = e L 12 T converge to the a value y * value as n increases for all AND gates and for two values of the temperature, and J 3 = 1 and J 12 = 0.1. The value, while being not one, is close to the phase transition value y * = 1, as we can see in Fig. 10. We interpret this as a symptom of metastability. This is the case for both horizontal and non-horizontal interactions. However, with horizontal interactions, the effective size of the system is larger due to the presence of the virtual spins. Thus, in the case with horizontal interaction the system reaches the ground state with more difficulty. For a larger number of layers (n ≥ 9), the zeros converge to a finite value for various values of the temperature. Another way of interpreting Fig. 10 is that in the case of homogeneous trees (all gates identical) one never effectively has a phase transition for a finite system, while it is the case for an infinite one; this is an important difference between the Cayley tree and the Bethe lattice [48]. While this is true for ferromagnetic trees (corresponding to AND for instance), it is not true for random gates [51], in which a second order phase transition is known to exist. Let us use for instance the case of the Kibble-Zurek scaling [49,50]: the density of defects scales as ρ = ξ −d , where d is the dimension of the system and ξ the correlation length at the freezing time, which is related to the linear annealing parameter τ as ξ ≈ √ τ . From the exponential annealing, we see that τ −1 = −T c log(λ)/J. For a quasi one-dimensional system, like the one of interest in this paper, we have thus Q = 1 − ρ 1− ≈ T c (1 − λ) as λ → 1 (i.e. infinitely slow annealing). Therefore, in our simulations in which we have random gates, we conjecture that the defects we observe are due to the transition observed in [51].

Conclusions
We proposed logical gates made of interacting spin and tested the possibility of integration in tree-like circuits. We find that in absence of interaction loops and for finite systems, it is possible to mimic deterministic computations by integrating logic gates made of interacting magnetic bits. In the case with loops, the presence of internal, faulty gates, best described in terms of temperature dependent probabilistic gates, suggests directions toward probabilistic [45] rather than deterministic computing. The computation of a specific Boolean configuration can be in general thought of as one particular configuration of random field on a tree, which presents interesting universal behavior [51] and will be studied elsewhere for our case. One of the drawbacks of our systems is that it scales exponentially with the number of layers, implying that only a of the partition function in the complex plane, from the recursion relation with horizontal interactions. We calculate the values of the zeros for J 3 = 1, and J 12 = 0.1, and we solve numerically for the location of the Yang-Lee zeros, as a function of the number of layers (n = 5 · · · 13). The maximum number of layers that we are able to resolve numerically is n = 13, which can be achieved thanks to the recursion relation. As the number of layers increases, the zeros move to the right and accumulate around y = 1. This would hint towards the emergence of a phase transition, also for small values of J 12 and faster for larger values. However, in Fig. 10, we see that the convergence of the zeros stops for a larger number of layers (exponentially many more spins) to a value which is not 1.  Figure 10. Convergence of the Lee-Yang zeros with the number of layers, and as a function of T for J 3 = 1 and J 12 = 0.1. We write z = e iθ , and θ → 0 would imply the emergence of a phase transition. We see that as the system increases (exponentially) in size, the zeros do not converge to the zero. In fact, as the system increases, the zeros converge to a finite number in the complex plane.
few layers might be practical. Another potential problem in practice is the local, biasing field, necessitating islands of different coercivity. Nonetheless, our analysis cautiously suggests the viability of deterministic computation using interactive, magnetic memory bits, realized by single domain magnetic nanoislands. In future work, we will explore the probabilistic computing aspects of these systems. In the equation above, we have introduced the scalar spin variables σ i and σ j . Let us consider now the case where all the spins are pointing upwards and focus on the case of the interaction between the bottom and the top spin of the triangles. In this case θ 1 = θ 2 . As a convention, the bottom spins are σ 1 and σ 2 , and the top spin is σ 3 . It is then easy to see that we can write θ 1 = θ 2 ≡ θ as a function of r and h, as θ = arctan( 2r l ), and thus: and as r → ∞, θ → 0. We see that there is a specific angular dependence on the coupling sign. We see that for r = r * = √ 2l 3 the interaction parameter is zeros, and for r > r * , the interaction changes sign. Regarding the horizontal spins, we have σ i · r 12 = 0 for i = 1, 2, and the energy (A. 6) In the bulk of the paper we exploit the fact that in the case of magnetic nanoislands, which possess a dipole moment, it is possible engineer the position and direction of the dipoles in order to choose the sign and the interaction strengths.

Appendix A.2. Gates with horizontal interactions
The first generalization we perform is to add horizontal interactions between spin σ 1 and σ 2 in the tree model introduced. As a first generalization, we consider the case in which vertical interactions and horizontal interactions have the same coupling. In the next sections we consider the case in which the spins interact on isosceles triangles on the plane. This means that vertical interactions are identical, but because of the property of the property of dipole interactions described above, the horizontal interaction can have a different value related to the relative angle between the spins. Such inclusion is also in view of plausible experimental implementations in which horizontal interactions are hard to suppress. We report for completeness the table of truths for the AND, OR, XOR, NAND and NOR in Tab. A1. The Ising model description of the interaction between the two models is given by: Table A2. We observe that the addition of the horizontal interaction implies a nontrivial shift to the energies.
which are presented in Tab. A2 as a function of the spins σ 1 and σ 2 . We are interested in the behavior of the spin σ 3 , which we interpret as the output of the interaction, given σ 1 and σ 2 . We write E(σ 3 |σ 1 , σ 2 ) as the energy of the Hamiltonian of eqn. (A.7) given the values of σ 1 , σ 2 , and as a function of σ 3 . We can thus construct the Table A3, which shows how the energy gap emerges. Let us assume that we interpret the interaction above as a gate where σ 1 , σ 2 are the inputs, and σ 3 is the result, where −1 is interpreted as FALSE and +1 as TRUE. We ask whether which type of gates can be encoded in the ground state.
We ask whether again, given the parameters J and h, it is possible to obtain such table of truth in the ground state. For instance, for the gate AND, we ask whether there is region in the plane (J, h) in which: and for the gate OR: and finally, the gate XOR: The result is the one discussed in the bulk of the paper, and presented in Fig. 4, while the energies are presented in Tab. A4. We observe that meanwhile the gate XOR conditions are not all feasible, the gates AND and OR are feasible. Both of them require J < 0, i.e. ferromagnetic interaction, and for h > 0 we obtain an AND gate, meanwhile for h < 0 we obtain an OR gate. The situation is instead inverted for J > 0. This suggests that the local field h can program the system to obtain a behave as OR or AND gate. This picture suggests the use of external field to reprogram, given the sign of J, the type of gate one aims to use. The approach we use works if we pin the input spins σ 1 and σ 2 . The key problem when all σ's are free to fluctuate with this approach in general, is that there is strong degeneracy, and the minimum state is given by −1, −1, −1 for J ≤ 0. We can necessary to try to split the degeneracy by considering two set of couplings as mentioned before J 12 for the coupling between spin σ 1 and σ 2 and J 3 between spins σ 1 and σ 2 with σ 3 . We thus consider the following Hamiltonian: and this Hamiltonian better represents the dipole interaction between not uniaxial dipoles. The energy states are presented in Tab. A4 forthe 8 states between the three spins.
In the previous calculation for the average of σ 3 , the output spin, we have used the couplings J's between the spins to calculate the average. The same can be done, however, by the mentioned star-triangle transformation between the spins, thus the couplings n's. As mentioned, this mapping introduces an unphysical spin σ which directly couples to the three physical spins, which we call now σ .
The effective Hamiltonian for the interactions is given by: where h is normalized in the temperature. Let us consider now the average σ 3 as a function of the temperature for these couplings. The average output spin is given by:

Appendix B. Partition function recursion
We are not interested in deriving a recursion relation for the partition function of the system both for J 12 = 0 and J 12 = 0, via the star-triangle transformation. The importance of such recursion relation is that it simplifies the numerical study of the Lee-Yang zeros. We use the formalism introduced for the star-triangle transformation of the physical σ and the unphysical σ spins. We define Z(σ r , h ) = σ r e h σ rZ (σ r ). We also consider for simplicity h homogeneous for the time being. If the root is the spin σ 0 , then we have: = 2 cosh(h +L 3 )Z(σ 1 ,h +L 12 )Z(σ 2 ,h +L 12 ) It is easy to see that such recursion relation can be generalized. It is simply necessary to replace h with h σ j . Thus we have: Figure B1. The recursive structure of eqn. (B.1): Z n (σ, h) represents the partition function rooted at the spin σ with external field h, which depends on the partition functions at the lowest order Z n−1 rooted at the two sub-branches spins.
If we assume thatL 12 = 0, then it further simplifies to: If we introduce the variable y = e h , then it is easy to see that cosh(h σ 0 +L 3 ) can be written as Z n (σ 0 , y) = (1 + y 2 ) cosh(L 3 ) y Z n−1 (σ 1 , y) 2 . (B.5) Let us assume that the zeros of Z n−1 are y n−1 1 , · · · , y n−1 k , and thus Z n−1 = a n−1 k i=1 (y − y n−1 i ). We have: and thus we have that the zeros at the order n will have two further zeros: which are imaginary. We note that this is independent from the sign ofL 3 . In the case the external field is not homogeneous the result easily generalizes. In fact, we have: and thus if Z n−1 does not have any real zeros neither will Z n . The result follows from noticing that Z 0 (y) = 1+y 2 2y , which does not have any real zeros. We thus see that the model with L 12 = 0 cannot have any phase transition.
Appendix B.2. CaseL 12 = 0, L 3 > 0 Let us now consider the more interesting case in which L 12 = 0. In this case, the model is not exactly a tree anymore.
From the combinatorial point of view, the terms which appear in the partition function can be obtained from the analysis of the tree in