Artificial Spin Ice Phase-Change Memory Resistors

We study the implications of the anisotropic magnetic resistance on permalloy nanowires, and in particular on the property of the resistance depending on the type of lattice. We discuss how the internal spin configuration of artificial spin ice nanowires can affect their effective resistive state, and which mechanisms can introduce a current-dependent effect dynamic resistive state. We discuss a spin-induced thermal phase-change mechanism, and an athermal domain-wall spin inversion. In both cases we observe memory behavior reminiscent of a memristor, with an I-V hysteretic pinched behavior.

Here we explore the possibility of using artificial spin ice to engineer memristors, based on previous work on transport and magnetoreristance in these materials [49][50][51], and start from the theoretical framework previously developed by one of us [52]. It was shown there that connected artificial spin ice is as a electrical circuit where tension drops at the vertices because of the magnetoresistive effects of domain walls there. Vertices can be considered electrical elements whose functionality is controlled by the magnetic moments impinging in them: changing their configuration affects the resistance of the system, leading to reconfigurable circuitry. However one can imagine that current itself could alter the moment configurations, thus leaving memory of its passage. We consider two possibilities for this coupling. One is based on Joule effect, where the superparamagnetic temperature of the nanoislands exceeds by little the operative temperature of the system [53]. The other is through the spin torque of a spin polarized current [54][55][56][57].
In this work we explore whether connected artificial spin ice can function as a memristor by solving the collective dynamics of currents that alter the magnetic texture, which in turns alters the localized resistance and then the currents themselves. We first derive a perturbative equation for the effective resistance of the device as a function of the nanoislands moment configurations. We then use this exact solution to obtain via self-consistency a closed equation for the conductance of the device. This latter can be written as the sum of the conductance of the permalloy nanowires and a state dependent conductance. We show first the simple case of a 3-legs junction, and then show that in general for materials with a sharp magnetic order transition the resistance is of the phasechange type. We then extend our study by simulating a Kagome lattice [58,59] when a threshold domain wall spin inversion is considered, and observe the effect of the many-body interaction on the resistance.
FIG. 1. The structure of the anisotropic magneto-resistance memristor device we study, within a Kagome lattice. We consider a system of resistive nanowires of magnetic permalloy, which are connected to a battery. The currents flow in the nanowire, but because of the anisotropic magnetoresistance induced by a small external magnetic field H, the internal resistivity depends on the distribution of the magnetic moments in the wires at the junctions.
Consider a network of permalloy nanowires, as shown in Fig. 1. Each wire portion between vertices is magnetized and the coupling between moments is such to obey the ice-rule. If the external magnetic field is zero, each wire of the network will have a resistance R = ρ 0 L, where L is the length of the island. The presence of magnetization m alters the resistance according to the Anisotropic MagnetoResistance (AMR) law at low magnetic field, H where J is the density of current, and ρ || (H) and ρ ⊥ (H) are the resistances parallel and perpendicular to the magnetization, respectively (and ∆ρ = ρ || −ρ ⊥ ). Along a line γ in the material, the voltage drop is and as will shall see it affects the behavior of the device by introducing a state-dependent memory. As noted previously by two of us [49,52], the AMR is independent from a global change of spin configuration. If V 0 is an external potential applied to the circuit in one direction (see Fig. 1), the current follows Ohm's law, where R m is the resistance of the galvanomagnetic material of interest and R v is an external resistence. We aim to show that the internal configuration of the magnetization acts as a memory for the resistance R m . We neglect any parasitic capacitance. We assume that the magnetization of the system is at a equilibrium temperature at a certain temperature T .
Two of us showed previously [49,52] that the voltage drop across nodes depends on the configuration of the moments in the nanowires {s i }. The problem then becomes how to obtain the current distribution when the magnetic moments configuration is known. For that we employ a graph theoretical approach [38,60], which has been already successful in the study of circuits of memristors [38,39,61,62].
Consider a general graph G (for definitness a Kagome spin ice in Fig. 1) with N v vertices (or nodes) and N e edges, which describes a network of resistors. The graph supports N c loops (closed loops or subcircuits). One can describe the potential equivalently by assigning a potential p α at each node, or a potential drop v k for each edge k hosting a current i k (we use latin indices for the edges, and greek indices for the nodes, greek indices with tildes represent instead cycles on the graph). We choose an orientation for the current on each edge-something that can be done in 2 Ne ways-and encode it into the matrix In this language, M j=1 B αj i j = (B · i) α = 0 enforces Kirchhoff current law at each vertex α. Then the potential drop v k for each edge k along the chosen direction is The Kirchhoff Voltage Law (KVL) can be written as k Aξ k v k = 0 where Aξ m is the N c × N e cycle or loop matrix, obtained by assigning first an orientation to the edges and then the loops of the graphs, and assigning values of +1, 0, −1 if the orientation of the loop agrees or disagrees with the orientation of the edge. This equation simply states that the circuitation of the voltage on voltage on a node must be zero. As a consequence, in general, B · t A = A · t B ≡ 0. Finally we introduce R the N e × N e diagonal matrix of resistances.
As shown in ref. [38,60] one obtain a generalization of the Ohm's law in the form where we defined the symmetric matrix Q = −A t (ARA t ) −1 A. The vector V is the vector of effective voltages in series to the resistances on the graph (see the Supplementary Material, SM-C, for details). Its evaluation can be carried out analytically given a certain lattice configuration. In the SM-D, an analytical expression in implicit form is provided for the Kagome lattice, which will be used shortly. But first we discuss the implications of eqn. (3) when an external voltage is applied to the material, as in Fig.  1. The internal currents depend on the voltage, which in turn depends on the spin configuration within the material, and linearly in the currents as from eqn. (2). Then, a self-consistent nonlinear equation for the voltages can be obtained. Since the anisotropic magnetoresistance is a small effect (it contributes 3-5% on the material resistance) we can linearize Eq. 3 in ∆ρ. Remarkably, the final equation for the effective resistance is simply a parallel resistor equation (see SM-C and SM-D): where G 0 is the conductance when no anisotropic magnetic resistance is present, and G m (s) is a statedependent function with the dimension of an inverse resistance. In particular, G m (s) = ∆ρ Q t M (s) Q is a quadratic form where Q is a network dependent vector with the dimension of inverse resistance and can be obtained from the matrix Q from the rows (or column) corresponding to the resistance in parallel to the generator and the internal resistances. Instead, M (s) is an adimensional matrix which depends on the internal magnetization state via terms of the form s i s j . It is hard to obtain exact expressions for Q, but this can be easily obtained numerically. Naturally M (s) changes in time with the currents, leading to a memory effect that we aim to show to be memristive. How it changes depends on an underlying physical mechanism. We posit that if the nanoislands are close enough to the superparamagnetic threshold they can become thermally active due to teh Joule effect of the applied current. For illustration we consider a simple 3-moments system, as in Fig. 2. In this case, the directionality of the moment now take a rather small set of possible resistances which introduces different resistive states. For larger system however we can introduce an effective description for the resistance.
We can write an approximate (first order contribution) equation for the thermal average of the effective resistance of the form: where R < and R > are two resistances which depend on the value of M (s) and the geometry of the spin ice (which is contained in Q). The resistances are defined as quadratic forms (see SM-C) where resistances can be obtained also if the interaction is antiferromagnetic (and the graph bipartite). The smoothed theta function, in this case, is a reasonable smoothed θ kfunction is θ k (x) = 1 1+e −kx , where θ(x) = lim k→∞ θ k (x) which incorporetates a sharp transition or a crossover. The effective mechanism is based on the fact that that the equilibrium temperature depends on the balance between radiation and current induced Joule heating. Thus, in the simplest possible approximation, we see that under the application of a small magnetic field the resistance of the material can be assumed to be a thermal phase-change type of material (or switches), in which the system has two resistance phases depending on the current, which in turn controls the temperature of the permalloy via Joule heating. If the transition is not sharp but only a crossover, then we can assume that whereR is a smooth function such thatR → R < as T → o and andR → R > as T → ∞. In this case, it seems reasonable to assume that the material will fall in the thermal memristor framework introduced in [53]. Then, because we expect the typical memristive v − i hysteresis to be small, it is possible to see the change in the resistances from the v − r Lissajous figures (Fig. 2), obtained from the functional dependence of the effective resistance which we have obtained.
Another mechanism for moment inversion in nanowires which does not require very careful fine tuning of the temperature, is the current-induced domain wall inversion via spin-transfer [54]. Specifically, when a current is applied to a magnetic nanowire, some domain wall defects can form at the junctions and quickly move along the wires [55][56][57]. This phenomenon can be effectively modeled in our system by assuming that if the current in a wire is higher than a certain threshold I c , then the spin along that wire is inverted.
In a large system this hard switching behavior gets smoothed out. We consider an extended Honeycomb circuit of 16 loops, where heach whire has resistence R 0 . We then consider the following two-steps spin dynamics. At each time step we start with a spin configuration s t−1 , solve for the Kirchhoff laws in the nanowires, and find the equilibrium currents i(t) in the material as a function of the external voltage V (t) = V 0 cos(ωt) only. We then use the eqns derived in the SM-C to find the auxiliary voltage generators (given the spin configuration s t ) in the material which will affect the spin configuration i(t). At this point, we consider the domain-wall inversion process. If the current in each wire is above a certain threshold I c and the current is in the opposite direction of the magnetization direction of that wire, we flip the magnetization direction instantaneously. In experiments like those of [54] the switching is extremely fast. We thus consider for simplicity instantaneous inversion.
We account for the manybody interaction among spins by imposing constraints on the possible vertex configurations. Honeycomb spin ice [52,58,59,63] is frustrated and at low energy enters an ice-rule regime where only vertices with two moments pointing in and one out, or viceversa, are allowed. Therefore we only allow spin inversion when it produces vertices of magnetic charge equal to Q = ±1, and neglect the nodes with Q = ±3.
Once we fixed all the parameters, there will be a threshold V c which depends on the size of the system and the resistivity of the material above which the spininversion occurs and below which the system is a normal resistor. The threshold dynamics is reminishent of a fusenetwork dynamics [64], but with the difference that it is not the conductance that dramatically drops, but that the effective voltages change instead.
In Fig. 4 we plot resultsof current vs. voltage. For small values of V 0 no memristive behavior is present (inset), as expected. For both the non-interacting (red) and interacting (blue) honeycomb lattice we obtain a zerocrossing pinched hysteresis loop typical of memristive devices and which suggests the presence of memory [41][42][43]. The latter is more hysteretic and smoother than the former but memristive behavior is already present even in absence of interactions. The area of hysteresis is small due to the small value of the magneto resistence effect.
We have put forward a theoretical framework for the study of memory properties of magnetic nanowires subject to an external field as an effect of anisotropic magnetoresistance, building on previous results [52]. We have derived exact and general equations which show, given a certain spin ice lattice, how the resistance of the material changes given the internal configuration. This has enabled us to obtain first order contributions to the resistance, showing that there exist an effective resistance in parallel to the nanowires network and which depends on the internal state of magnetization. Then teh coupling between current and magnetism lead to a memristive behavior. We studied two mechanisms which induce these memory effects, and which likely coexist. As the anisotropic magnetoresistance effect is of the order of a few percentages, we expect a smaller change in the value of the resistance as a function of the voltage. However, we have shown that the hysteresis curves depend on the many-body interaction and the configurational manifold of artificial spin ices. This suggests that more generally the functionality of artificial spin ice memristors can be open to design-as so many other properties of these materials proved to me. These ideas are an alternative to Spin-Torque memristors for bio-inspired computing [65], to produce an effective magnetic phase-change material [66]. In fact, because of the sensitivity to temperature, memory resistors in spin ices can have a variety of behaviors that can serve as an alternative to known Spin-Transfer-Torque devices [65,67] and Phase-Change Materials [66]. Furthermore, as the magnetic moments can be acted upon collectively by external magnetic field, or individually, artificial spin ice memristors could be reprogrammable.
The work of FC and CN was carried out under the auspices of the US Department of Energy through the Los Alamos National Laboratory, operated by Triad National Security, LLC (Contract No. 892333218NCA000001). CN was founded by DOE-LDRD grant 2017014ER. FC was also financed via DOE-LDRD grants PRD20170660 and PRD20190195.

A. Formal solution of linear circuits
We use a graph theoretical approach [1,2] to solve for the current knowing the nanoisland moments. Consider a graph G (for definitness a Kagome spin ice in the figure) with N v vertices (or nodes) and N e edges, which describes a network of resistors. The graph supports N c cycles, that is closed loops or subcircuits. In each node there is a potential p α , and for each edge a current i k (we use latin indices for the edges, and greek indices for the nodes, greek indices with tildes represent instead cycles on the graph). We choose an orientation for the current on each edge-something that can be done in 2 Ne ways and encod it into the matrix B αk of size N × M . Then M j=1 B αj i j = (B · i) α = 0 enforces Kirchhoff current law at each vertex α. Then the potential drop v k for each edge k along the chosen direction is given by v k = ξ p ξ B t ξk = ( t p · B) k . The Kirchhoff Voltage Law (KVL) be written as k Aξ k v k = 0 where Aξ m is the N c × N e cycle or loop matrix. This equation states that the circuitation of the voltage on voltage on a node must be zero. As a consequence, in general, B · t A = A · t B ≡ 0. In order to see how formally one can introduce the reduced loop matrix, such that (ARA t ) is invertible, we need some notion of graph theory. Given the graph G, we introduce a spanning tree T (called co-chords), and the set of edges of the graph not included in the tree as T , or chords, are given byT . For each element of the chordT , we assign a cycle, called fundamental loop. The loop matrix A, can be reduced to its N e − N v + 1 fundamental loops. Then, it is not hard to show that the current vector can be written as i = (A t T i c , i c ) = A t i c , where we used the fact that given a chords and co-chords splitting, we have ( Since A is derived from the reduced incidence matrix, this is called reduced loop matrix. At this point, it can be shown that which is the starting point of the paper. It is not hard to see that ARA t is always invertible for non-zero resistances. For more details, we refer to [1,2]. The reduced loop matrix A is constructed using the following procedure. First, we assign an orientation to the edges of the graphs, and for each loop of the graph, we assign an arbitrary orientation to the loop along each edge of the loop. We then construct the matrix A LE (dimensions of loops by edges) as follow. If the loop N c does not contain the edge E, A LE = 0. If the orientation of the loop agrees with the orientation of the edge, then A LE = 1, otherwise A LE = −1. At this point, we choose a subset of N e − N v + 1 linearly independent loops and remove the others from A. What we obtain is the reduced loop matrix.

B. Mapping voltages drops at nodes to voltage generators
It is common in spin ice materials to approximate the magnetization with a internal configuration M = { s i }, where s i = s i {ax, bŷ} are Ising variables on the plane, which cannot rotate. The plan is to map the node configuration to a set of voltages in series to the resistances, as this is an exactly solvable model For each edge β (which represents a resistance) between the nodes (n i , n j ) and, and we consider a tuple of voltages (V β,i , V β,j ) associated to it as in Fig.  1. Let us call F i the number of resistances attached to the node n i , which in graph theory is commonly called degree. We define also V β,i The goal of this section is to derive the voltages V β,z based on the configuration of the spins, which as we will see is connected to the voltages E i β1,β2 below. As introduced in [3], the node configuration can be assessed via the voltage integral across the node, starting from a resistance β 1 and going into a resistance β 2 . In the formalism of the anisotropic magnetoresistance, given a certain local node n i , and spin configuration at that node, a number of voltages E i β1,β2 can be obtained via the anisotropic magnetic effect: Let us call G the graph that represents the circuit.

Bulk
If the graph is planar, then if K is the number of resistances entering a node i, because of the planarity of the graph, only K values of E i β1,β2 for a given node. This is due to the fact that for planar graphs only a number of cycles equal to the number of faces of the dual graph are necessary to obtain a self-consistent equation. However, the number of faces in this cases equals the number of entering edges. Thus, a very natural choice is to choose a set of fundamental loops in the circuit that are associated to each node in the dual graph. Also, because the graph is planar, we can choose a consistent orientation for each (fundamental) cycle in the circuit. Given this prescription for each node n i , the number of integrands E i z,z+1 is equal to the number of voltages V β,i . In particular, we have the relationship, obtained by performing the integration via eqn. (2), and the voltages V β,i . For each node n i , let us call B i the set of edges incident to that node. Because of the planarity, it is possible to give a consistent ordering to the edges B i = {b i 1 , · · · , b i Fi } as well such that b r+1 − b r = 1, as in Fig. 2. Then, we have where F i and F j are the number of resistances attached to the nodes n i and n j respectively. It is not hard to see that the equation above, for each node, can be written in the more compact form: where the matrix F D is a matrix of size F × F given by: which is clear to be the discrete derivative on a circle with F points. Thus, it is clear that for each node, this matrix is not invertible as it contains one null eigenvalue, with eigenvector proportional to e i = (1, · · · , 1) t of arbitrary sizes for an arbitrary constant c i associated with each node, and where we calledD −1 i the pseudo-inverse operator. As we will see however the choice of this constant does not have any physical implication and we can freely set it to zero. The pseudo-inverse for the forward difference operator can be written as FD −1 = ( FD t FD ) −1 F D, where FD −2 is the pseudo-inverse of the second difference operator. We focus for now on the pseudo-inverse FiD −1 of the matrix, that can however be explicitly calculated, as we know that F D F D t = F D t F D = F D 2 , e.g. the discrete second derivative on the circle of dimension F . Thus, for each spin configuration of the spin ice at each node, given by the associated voltages E i , we have a vector of the effective voltages on each edge β, which can be written as where q β = c i − c j .

Fixing of boundary resistances
The inversion problem at the boundary is slightly more complicated than the one in the bulk, and requires the prescription of setting some voltages to zero to avoid overdetermination. Given a circuit with a well-identifiable boundary, the it is not hard to see that given the prescription of Fig. 1, the system of equations is underdetermined. At each node on the boundary with m resistances, we have m − 1 loop constraints. Thus, we need to find a way to reduce the number of voltages on the boundary for each node. Our prescription is the following. Let us consider the set of boundary resistances R b = ∂G. Because the graph is assumed to be planar, we can assign a consistent orientation to this boundary, O. Then, our prescription is that, given V β,k if the orientation of the boundary and the positive side of the voltage generator agrees, then we keep it, while otherwise we remove it (or set it to zero). For instance, in Fig. 3, given the orientation of the boundary (red arrow), the generator highlighted in red is set to zero. It is not too hard to see that this prescription removes the extra degree of freedom at each node on the boundary.

C. General approach: absorbing the spin configurations in voltage sources
As it is shown in the Appendix, given the node dependent voltage configurations of eqn. (eq:amr), we can obtain equivalent voltage generators depending on the configurations of the spins V ({s i }). This is important, as we can now write an exact equation for the currents of the system at equilibrium, as this is a resistive system with voltages in series. The solution is known and given by: where A is the directed cycle matrix on the fundamental cycles of the circuit. Here we assume that the voltage V is indeed depending on the internal spin configuration s.
We now comment on the constants c's. It is interesting to note that these can be written as q = B t c, where B t is the directed incidence matrix of the graph. However, it is known that AB t = 0, and thus any configuration of these constants has no impact on the configuration of the currents, as one would expect from a change in potential. Another way to see this is by noticing that for each fundamental loop, necessarily at each node the same constant must be counted twice. However, since the cycle is directed, via the Kirchhoff law the same constant appears twice but with opposite signs, as it can be seen in Fig. 2. We can thus set these to zero.
Let us now discuss how to write the effective memory of the component. The voltage V is n + 1 dimensional, where n is the number of edges internal to the device, and 1 is the edge where the external voltage is applied. First, let us call the matrix Q = −A t (ARA t ) −1 A. The diagonal matrix R can be written as diag(r, · · · , r, R v ), where r is the resistance of a single alloy nanowire, while R v is the resistance of the battery. For the matrix Q, we consider the following splitting: which is necessary to distinguish the resistance of the device, and the resistance of the battery. Let us call ( V ) 0 = v 0 the applied voltage, and the rest n-dimensional vector V r , which are the internal voltages. Similarly, we introduce the splitting of the currents i, as ( i) 0 = i 0 and i r , the n−dimensional vector associated with the internal currents. Clearly, at equilibrium, we must have that these voltages depend linearly on the magnetic anisotropic effect, and on the internal configuration of the spins. We can write, because of eqn. (2), The equation above can be written, given the splitting of eqn. (9), as The internal currents at equilibrium are thus: Using the equation above for internal currents, we have Thus, at the first order in ∆ρ, we obtain that It is not hard that we can re-write eqn. (14) in terms of resistances. We have where R 0 is the resistance when no anisotropic magnetic resistance is present. It is also interesting to note that the contribution to the conductance is a quadratic form. We see that the formula above states that the effective conductances due to the magnetic anisotropy and conductance of the alloys at zero external field sum. We note that however M (s) can change in time due to the currents. Thus, the equation above states how the effective resistance changes with the internal degrees of freedom. For h → 0, the effective resistance due to the magnetic anisotropy goes to infinity, and since these are in parallel, the resistance of the material goes to its original value. It is thus now the goal to construct the matrix M (s).

D. A worked out example: mapping of the spin configuration to the voltages of Kagome ice
Let us consider the case of nodes with 3 legs, in the case of the Kagome lattice in Fig. 5.  As we have seen, this can be done node by node. For each node, we are going to have, given a configuration of the three spins incoming to that node, a voltage configuration which depends on the spins s 1 , s 2 , s 3 . For clarity, s i is positive if it points towards the right We know that the voltage configuration is independent from a change of the sign of the three spins, as this results in a change of direction of the magnetization at the node, and the voltage drop is independent with respect tom → −m. For each edge, we are going to have four possible configurations of the spins at the node: Thus, given a 3-dimensional vector E for a node for each of the four configurations above as: where δ s is a Kronecker delta which is one if s = 1 and zero if s = −1. This Kronecker delta can be written δ s = 1−s 2 .
It is now not hard to see that we have: On a Kagome lattice we have two type of nodes. Let us call them 1 → 2 and 2 → 1, as in Fig. 6. We assume that in both cases the currents direction are from the left to the right, that the integration over the cycles are clockwise and that the spins are positive if they point right. We work first with the 2-1 node. We call s 1 and s 2 the in-nodes and s 3 the out-node.
In this case, FIG. 6. The 2-1 and 1-2 nodes, and the associated direction of the node magnetization for each associated spin configuration.
In the case 1 − 2 instead, we have The factors of 1 2 come from the projection onto the current directions (cos( π 3 ) = 1 2 ). The effective magnetic moment at the node is in fact either directed towards the link, or has an angle π 3 . We stress that the signs of the currents depend only on the direction of the integration of the voltage over the node with respect to the direction of the currents, and not on the magnetization.
The voltage at each link, since it is the difference of two voltage sources, depends on five different spins, which decide the magnetization of the nearby nodes.
We have that i From the expression above we see that M (s) is not-diagonal, but that for the Kagome lattice is a block which involves 5 currents. One immediate example is an horizontal edge in the ground state. The non-zero portion of the matrix M (s) which corresponds to the voltage V 3 and the currents i 1 , · · · , i 5 as in Fig. 7 is given by: We see from the expressions above that this formalism is ought to be used for a numerical simulation rather than for analytical computations, and that M (s) depends on pairs of variables s i s j .

Magneto-resistance memory
We now consider the simplest non-trivial example of magnetoresistance memristor device. In this simple example we focus on a simple enough case for which most of the techniques we have developed for general constructions are not necessary. In particular, only one voltage per edge is necessary, and thus simply the extra voltage vector can be simply written as and if we use the equilibrium current equation: from which we obtain Upon investigation, we find that the product of the matrices M (s) andD −1 are given by: which we will now use. We see that the matrix which couples the internal spins to the internal currents is a rather non-trivial matrix which, however, depends only on the internal configuration. In the next section we show that when the spins are allowed to flip thermally, a non-trivial memory effect ar ises out of equilibrium. Thus, the effective voltage in this case is simply As a result, we have which is the state-dependent effective voltage.

E. Thermally induced flips: out of equilibrium properties
It is interesting at this point to observe the out-of-equilibrium dynamics of the system. We perform numerical simulations, apt at enhancing the effect and to show how a hysteresis loop typical of a memristive system emerges in this scenario.
At equilibrium, the currents satisfy the Kirchoff laws. For ∆ρ = 0, the system does not present any difference from a normal resistance. However, as ∆ρ = 0, thermal coupling can affect the internal properties of the resistance.
Here we assume a very simple internal dynamics, governed by the thermal coupling due to the Joule heating of the device. The model we suggest is rather simple but explicative of the phenomenology. The internal state of the device is assumed to evolve according to a Metropolis dynamics for the 3 spins, with a flipping probability: where ∆H is the energy difference between one configuration and the proposed one. The energy is the simplest possible ferromagnetic coupling for 3 nanoislands, given by The coupling between the internal states and the currents occurs via Joule heating: as the currents flows, we assume a temperature for the devices which follows a very simple relationship: where the first term is due to the Joule heating effect (and we assume that the temperature is just the average heating of the three branches), and as a balancing effect for the temperature we consider Stefan's radiation law. Given this simple mechanism, we consider the out-of-equilibrium voltage v(t) = v 0 sin(ωt), which is shown in Fig.  Fig. 8. We observe zero-crossing hysteretic jumps due to the state of the spins, between two resistance lines. In order to observe a real memristive behaviour, we need to go to larger systems.

Effective model after thermal averaging
For a larger lattice, obtaining the matrix M (s) can be challenging. However, the key features of the resistance can be inferred from thermal averaging as follows. Let us consider eqn. (14) again: In particular, we are interested in the thermal average of the equation above, e.g.
where · T is the thermal average, over all the possible configurations of the system at temperature T . We note that is a matrix which, for a local Hamiltonian, is composed of products of neighboring spins only, of the form where the distance d(s i , s j ) is of order one. Depending on the system of interest, this average will lead to different results depending on the geometric arrangements of the nanoislands. In particular, for lattices which exhibit with a sharp transition from a disordered to an ordered phase at low temperature, we can approximate Thus, because the matrix M (s) is composed only of products of pairs of neighboring spins, we can write and we can think of an effective interpolation between two limiting values of the resistance. Given this feature, can write an approximate (first order contribution) equation for the thermal average of the effective resistance of the form: where R < and R > are two resistances which depend on the value of M (s) and the geometry of the spin ice (which is contained in Q), and are defined by The function θ k (x) is a smoothed Heaviside-theta function. Thus, in the simplest possible approximation, we see that under the application of a small magnetic field, and of joint permalloy islands (nanowires), the resistance of the material can be assumed to be a phase-change type of material (or switches), in which the system has two resistance phases depending on the current, which as a matter of fact controls the temperature of the permalloy via Joule heating. If the transition is not sharp but only a crossover, then the it is not too daring to assume that whereR(T = 0) = R < andR(T = ∞) = R > is a smooth function. In this case, it seems reasonable to assume that the material will fall in the thermal memristor framework introduced in [? ]. In this case, because we expect the typical memristive V − I hysteresis to be small, it is possible to see the change in the resistances from the v − r Lissajous figures, obtained from the functional dependence of the effective resistance which we have obtained. Albeit the exact numbers will depend on the type of material, we expect to be able to distinguish the type of transition from the V − R curves of the device as a function of the frequencies. At slower frequencies, the changes in the resistance will be more symmetric in the continuous case, while more abrupt but still hysterestic in the discontinous case. * caravelli@lanl.gov † gc6u@virginia.edu ‡ cristiano@lanl.gov