Electrically programmable magnetic coupling in an Ising network exploiting solid-state ionic gating

Two-dimensional arrays of magnetically coupled nanomagnets provide a mesoscopic platform for exploring collective phenomena as well as realizing a broad range of spintronic devices. In particular, the magnetic coupling plays a critical role in determining the nature of the cooperative behavior and providing new functionalities in nanomagnet-based devices. Here, we create coupled Ising-like nanomagnets in which the coupling between adjacent nanomagnetic regions can be reversibly converted between parallel and antiparallel through solid-state ionic gating. This is achieved with the voltage-control of the magnetic anisotropy in a nanosized region where the symmetric exchange interaction favors parallel alignment and the antisymmetric exchange interaction, namely the Dzyaloshinskii-Moriya interaction, favors antiparallel alignment of the nanomagnet magnetizations. Applying this concept to a two-dimensional lattice, we demonstrate a voltage-controlled phase transition in artificial spin ices. Furthermore, we achieve an addressable control of the individual couplings and realize an electrically programmable Ising network, which opens up new avenues to design nanomagnet-based logic devices and neuromorphic computers.


Main text
The ability to electrically manipulate magnetism is crucial for spintronic applications including spin-based data storage and computation devices.The successful switching of the magnetization in nanomagnets by means of spin transfer torques 1 , spin-orbit torques [2][3][4] , voltagecontrolled magnetic anisotropy (VCMA) [5][6][7][8] and magnetoelectric coupling 9,10 have led to significant steps towards next-generation low-power and high-speed magnetic memories.However, while nanomagnet-based data storage relies on this magnetization switching, nanomagnet logic requires the engineering of the magnetic coupling that aligns the magnetization of adjacent nanomagnets with a specific relative orientation [11][12][13][14] .In addition, with the possibility to construct various two-dimensional (2D) coupled nanomagnetic networks, often referred to as artificial spin ices 15,16 , the control of lateral coupling is of particular interest for the investigation of collective phenomena such as magnetic frustration, emergent magnetic monopoles and phase transitions [15][16][17][18][19][20] , as well as for the implementation of multiple computation tasks such as Boolean logic operations [11][12][13][14]21 and neuromorphic computing [22][23][24][25][26] . Desite the different lateral coupling mechanisms available, including long-range dipolar coupling and nearest-neighbour chiral coupling, these are engineered either through the geometric design 27- 30 or by locally tuning the magnetic anisotropy during the fabrication process 19 .As a consequence, the functionality of coupled nanomagnets is determined once fabricated and the subsequent electrical modulation at run-time remains challenging.
Here, we present a mechanism to realize electrically tunable lateral coupling between two adjacent Ising-like nanomagnets by exploiting the VCMA-mediated competition of the symmetric exchange interaction and Dzyaloshinskii-Moriya interaction (DMI) that favour the parallel (P) and antiparallel (AP) alignment of the nanomagnet magnetizations respectively (Fig. 1a).Employing this concept for extended networks, we are able to explore the voltagecontrolled phase transition in artificial spin ices and to obtain an electrically programmable nanomagnetic Ising network that can serve as a neuromorphic computing element.

Voltage-controlled magnetic coupling
The voltage-controlled nanomagnets are fabricated from a Pt/Co-based multilayer that has a large interfacial DMI 31,32 and tunable perpendicular magnetic anisotropy 33 .As shown in Fig. 1b, the structure comprises two protected regions (in red) with a fixed out-of-plane (OOP) magnetic anisotropy and a 50 nm-wide gated region (in blue) with tunable magnetic anisotropy.
In the gated region, the Co layer is exposed to the electrolyte GdOx, constituting a solid-state ionic gate structure [34][35][36] and thus its magnetic anisotropy can be reversibly modulated between in-plane (IP) and OOP by applying positive and negative voltages to the top gate electrode, respectively.In contrast, in the protected region where an 8 nm-thick SiO2 layer is inserted between the Co layer and GdOx, the migration of oxygen ions from GdOx to the top interface of Co is blocked and hence its magnetic anisotropy is protected from voltage modulation.The detailed fabrication process is described in Methods and Supplementary Information S1.
The GdOx was grown in oxygen-deficient conditions such that, in the as-fabricated device, oxygen ions at the top interface of the Co layer in the gated region are partially absorbed by GdOx.The gated region exhibits IP magnetic anisotropy whereas the surrounding protected regions have OOP magnetic anisotropy, giving an OOP-IP-OOP anisotropy configuration.
Moreover, as a result of the interfacial DMI at the Pt/Co interface, the magnetization within the OOP-IP and IP-OOP transition regions twists following a left-handed chirality enforced via chiral coupling 19 .As a result, the two OOP magnetizations in the protected regions are effectively AP coupled (Fig. 1a).In order to obtain the AP ground state, an oscillating and decaying magnetic field is applied perpendicularly to the devices, serving as the demagnetization protocol (Supplementary Information S2).As shown in Fig. 1c (left panel), for an array of coupled elements (Fig. 1b), the magnetization in the two coupled OOP regions exhibits either ↑↓ or ↓↑ alignment.After applying a negative gate voltage VG of -2.5 V for 90 min, the gated region exhibits OOP anisotropy thanks to the formation of Co-O bonds promoting perpendicular magnetic anisotropy 33 .Due to the collinear alignment in the OOP-OOP-OOP configuration, the DMI influence is reduced and the symmetric exchange interaction results in ↑↑ and ↓↓ low energy states.In other words, the two OOP magnetizations in the protected regions are P coupled.The conversion from AP to P coupling is demonstrated with the demagnetized configurations given in Fig. 1c.This coupling can be described as an effective exchange energy ( ) where S1 and S2 represent the direction of the adjacent Ising macrospins, which can only point ↑ or ↓, and J is the coupling strength that can be tuned between P (J > 0) and AP (J < 0).
In order to systematically study the voltage-controlled coupling, nanomagnet elements were fabricated on a Hall bar for the electrical detection of the OOP magnetization via the anomalous Hall effect (Fig. 2a).Due to the large difference in size between the protected and the gate regions, the Hall resistance mainly reflects the state of the protected regions.Following the demagnetization protocol, we recorded hysteresis loops, which consistently show the switching between three magnetization levels, corresponding to ↑↑ (↓↓) at large positive (negative) fields and ↑↓ or ↓↑ at intermediate fields (Fig. 2b).Note that the Hall resistance starts at the middle resistance level corresponding to ↑↓ or ↓↑, indicating that the demagnetized magnetic configuration is AP.We then record minor loops to quantify the coupling strength between two protected regions.The minor loops are shifted horizontally by the exchange bias field Hbias = 127±16 Oe (-130±16 Oe) when the loop starts from positive (negative) saturation fields, confirming the presence of the AP coupling (red loops in Fig. 2c).The AP coupling strength can be estimated from J = HbiasMVOOP = -2.5±0.3 eV, where M and VOOP are the magnetization and volume of the magnetic material in the protected region, respectively.After applying a negative gate voltage, the magnitude of Hbias gradually decreases and eventually its sign is inverted, indicating the conversion from AP to P coupling (orange and green loops in Fig. 2c).The magnitude of exchange bias with P coupling is similar to that of AP coupling, and the strength is estimated to be 2.2±0.3 eV.We then verified that the demagnetized magnetic configuration is P as the Hall resistance started from either the highest or the lowest resistance levels corresponding to ↑↑ or ↓↓ (Fig. 2c), confirming the AP-to-P coupling conversion.
Additionally, the magnetic coupling can be reversibly converted between AP and P by changing the VG polarity.As shown in Fig. 2d, we altered the VG polarity back and forth 13 times and Hbias changed signs accordingly.Note that the time for the first conversion (~90 min) is longer than the subsequent conversion times (~30 min), which is also reflected by the reversal of the blue and purple loops in Fig 2c .This could be related to the additional energy cost of the first detachment of the ions from their original positions as well as to the formation of ionic conduction paths that provide subsequent faster conversion 37 .
When the coupling J is close to zero, the magnetic configuration obtained after different demagnetization cycles exhibits a stochastic behaviour.As shown in Fig. 2e, we start from the P state where the percentage of obtaining AP alignment, AP, is close to 0. When a gate voltage of VG = 2.5 V is applied, AP gradually increases over time and approaches 1.The percentage of AP alignment follows the Boltzmann law and is given by: where kB is the Boltzmann constant and Teff is the effective temperature resulting from the demagnetization protocol (Fig. S3).AP is then regulated by VG that modifies the coupling strength, giving AP > 0.5 for J < 0 and AP < 0.5 for J > 0. We note that AP remains almost unchanged after disconnecting the device for 8 hours, indicating that the coupling is nonvolatile thanks to the ionic gating effect 8 .We then apply a VG alternating between ±3 V and observe a synaptic plasticity of AP that varies according to VG (Fig. 2e).We repeat the demagnetization process 130 times at J < 0, J ≈ 0 and J > 0 (at I, II and III in Fig. 2e), and the magnetic configuration alternates between AP and P randomly with different percentages (Fig. 2f), which enables a modulation of the correlation between neighbour macrospins as required for Ising-type probabilistic computing [38][39][40][41] .

Mechanism for the AP/P coupling conversion
To provide deeper insight into the mechanism of the AP/P coupling conversion, we consider a macrospin model and calculate the coupling strength in response to the magnetic anisotropy of the gated region.As shown in Fig. 3a, S1 and S2 represent the direction of the magnetization in the protected regions, while Sg represents the direction of the magnetization in the gated region.The total energy of this system is given by the sum of the symmetric exchange energy (Eex), antisymmetric exchange energy (EDM) and the magnetic anisotropy energy (Ean) of the gated region, which can be written as: where Jex and Deff denote the effective exchange energy and the DMI vector.Jex > 0 and Deff < 0 in Pt/Co with left-handed chirality.<i,j> represents all the possible combinations of the nearest-neighbour pairs of macrospins.Vg is the volume of the magnetic material in the gated region and Kg is the effective anisotropy constant of the gated region, which is experimentally tuned between IP (Kg < 0) and OOP (Kg > 0) by applying a gate voltage.The energies for the AP and P configurations are then: eff g g g g eff g g

2
; when ; when and 2 ex g g ex g g P ex g g g g ex ; when 2 ; when The S1-S2 coupling strength is then determined from the difference between the energies EAP and EP and is given by: J = (EAP -EP)/2.On increasing Kg, J increases from negative to positive, reflecting the AP-to-P coupling conversion (Fig. S5).Notably, if the gated region is strongly IP (Kg << 0), J ≈ Deff <0, whereas if the gated region is strongly OOP (Kg >> 0), J ≈ Jex > 0, so providing an intuitive picture for the AP/P coupling conversion resulting from the VCMAmediated competition between the symmetric and antisymmetric exchange interaction.
In order to quantify the coupling strength, we carried out micromagnetic simulations using the MuMax3 code 42 .The initial magnetization in the protected regions is set to be ↑↓ (or ↑↑) with the cell magnetizations in the gated region pointing in random directions.The system is then allowed to relaxed until a stable magnetic state is reached.In the ↑↓ configuration, we find that a Néel domain wall separates the two protected regions, whose width decreases with increasing Kg.In the ↑↑ configuration, two 90° domain walls form, which eventually vanish as Kg is increased (Fig. S7).The total energy as a function of Kg is shown in Fig. 3b.The main features of the EAP (EP)-Kg and J-Kg curves agree with those of the macrospin model (Fig. S5).
When Kg < 0, EAP < EP and J saturates at -3.0 eV.When Kg > 0, EAP > EP and J increases with g K (Supplementary Information S3).The experimentally obtained coupling strength (-2.5 eV for AP coupling and 2.2 eV for P coupling) are also in good agreement with the values estimated from the micromagnetic model.

Voltage-controlled phase transition in Ising artificial spin ices
As the magnetization in the protected regions can only point ↑ or ↓, it is possible to create an array of nanomagnet elements that mimic artificial Ising systems and display magnetic phase transitions that occur on changing the coupling strength.We first construct a one-dimensional (1D) Ising chain of nanomagnets by repeating an alternating structure of protected and gated regions in a line (Fig. 4a).The AP coupling in the as-fabricated Ising chain leads to antiferromagnetic (AFM) order on demagnetization (Fig. 4b upper panel).After applying VG = -2.5 V to the chain for 90 min, all of the couplings are converted to P and the chain of nanomagnets exhibits ferromagnetic (FM) ordering on demagnetization (Fig. 4b lower panel).
We also vary the width of the gated region wg and find that the degree of both AFM and FM orders gradually decreases with increasing wg (Fig. 4c).Here the degree of the magnetic order is evaluated by determining the nearest-neighbour correlation function: where N represents the number of the nearest pairs of macrospins and the sum runs over all nearest-neighbour pairs.<SiSi+1> is equal to 1 and -1 for perfect AFM and FM order in the 1D Ising chain, respectively.
We then construct a 2D Ising artificial spin ice (Fig. 4d).The demagnetized square lattice exhibits an AFM checkerboard pattern in the as-fabricated state (Fig. 4e) 19 .A phase transition from AFM to FM order can then be achieved by applying a negative gate voltage (Fig. 4e).
This transition can again be quantified by calculating the nearest-neighbour correlation function <SiSi+1>.By altering the VG polarity, the magnetic order is switched between AFM and FM with <SiSi+1> varying between positive and negative values (Fig. 4f and Fig. S9).The <SiSi+1> for FM order is smaller than that for AFM order, implying that the P coupling is weaker than the AP coupling.This could be due to the fact that the dipolar interaction becomes considerable in extended lattices and inhibits the formation of the FM order (Supplementary Information S7).
By extending this scheme to a more complex 2D network, it is feasible to construct a programmable Ising network whose couplings can be electrically adjusted between AP and P.
Such a network has applications for both Boolean and non-Boolean computing.In the Supplementary Information S11, we describe the realization of reconfigurable Boolean logic gates such as a controlled-NOT and a controlled-Majority gate (Fig. S14), in which a dual logic functionality can be implemented by controlling the polarity of the control voltage.Moreover, our approach offers an efficient way to build an Ising-type neural network whose vertices are magnetically coupled to each other and each coupling is electrically adjustable.Many combinatorial optimization problems that are ubiquitous in fields such as artificial intelligence, bioinformatics, drug discovery, cryptography, logistics and route planning, can be mapped to an Ising network with specific Hamiltonians 45,46 .The solution of such problems can be obtained by finding the spin alignment corresponding to the ground state of the Ising network.In spintronic-based neuromorphic computing schemes, however, the couplings between vertices in an Ising-type neural network are usually achieved with additional CMOS circuits or by resistive crossbar arrays 40,43,44 .To illustrate the capability of solving combinatorial optimization problems using a magnetically-coupled network, we experimentally solve a benchmarking Max-Cut problem in an eleven-spin Ising network (Fig. 6).
The Max-Cut problem is frequently used for circuit design and machine learning 47,48 , and is one of the most basic combinatorial optimization problems.In a typical Max-Cut problem, one starts with a system (a graph) in which a certain number of elements (the vertices of a graph) are related to each other by pairwise connections (the edges of a graph) with assigned weights.Therefore, the solution of Max-Cut problems can be obtained by relaxing a physical system to its ground state within a finite time.Different Max-Cut geometries can be implemented using our approach (Figure S12).However, due to the geometric limitations of a 2D nanomagnetic network and the short-range nature of chiral coupling, only nearest-neighbour connections are possible.More general network structures can be implemented by exploiting the long-range dipolar interaction 23,49 , or by electrically coupling distant elements using the spin-transfer torque effect 40 (Fig. S13).Moreover, the ability to electrically program the magnetic coupling permits the adjustment of the Max-Cut problem at run-time, so enabling hardware-level programmability of the solver.

Conclusions
We have shown that we can electrically tune the magnitude and sign of the lateral coupling between nanomagnets by taking advantage of the antisymmetric exchange interaction and modifying the magnetic anisotropy in a gate region between the nanomagnets.The change in magnetic anisotropy is a consequence of the electrochemical reaction localized at the interface and ion migration in the gate dielectric under an electric field.The time required for the modification is limited by the reaction rate and ionic mobility, and can be reduced, in principle, using a high-mobility ion conductor 50 or the electronic-version of the VCMA effect, which is capable of operating at GHz frequencies 7 .
Our approach offers the possibility to investigate collective phenomena, such as the coupling-dependent phase diagram 51 and phase boundaries of mixed FM/AFM Ising-like artificial spin ices 52,53 , as well as the exotic magnetic phase of a spin glass 18,20 .Moreover, we have provided proof-of-concept demonstrations of reprogrammable nanomagnetic Boolean logic gates and combinatorial Ising solvers, which will inspire future research on unconventional computing devices based on nanomagnets.

Device fabrication
Films of Ta (5 nm)/Pt (5 nm)/Co (1.5 nm)/Al (2 nm) were deposited on a 200-nm-thick SiNx layer on a silicon substrate using d.c.magnetron sputtering at a base pressure of <2.7×10 −6 Pa and an Ar pressure of 0.4 Pa during deposition.The Al layer was oxidized to induce perpendicular magnetic anisotropy in the Co layer using a low-power (30 W) oxygen plasma at an oxygen pressure of 1.3 Pa.The fabrication of the voltage-controlled coupled nanomagnets was carried out with electron-beam lithography.Continuous magnetic films were milled into the shape of the bottom electrodes with Ar ions through a negative resist (ma-N2401) mask.
The upper Co/AlOx layers were milled through a positive resist [poly(methyl methacrylate), PMMA] mask to create the nanomagnets and lattice structures.Using electron-beam evaporation, a protective layer of Cr (2 nm)/SiO2 (8 nm) was deposited, with the protected region defined using a lift-off process through a second PMMA mask patterned by electronbeam lithography.Then an electrolyte layer of GdOx (30 nm) was deposited using reactive magnetron sputtering at an Ar pressure of 0.4 Pa and with a mixed gas flow of 50 sccm Ar and 1 sccm O2.In order to promote the ionic gating effect, a short milling process was implemented to partially remove the AlOx layer in the gated region.Finally, top electrodes of Cr (2 nm)/Au (3 nm) were fabricated using electron-beam lithography combining electron-beam evaporation with a lift-off process.The base pressure for the electron-beam evaporation was <1.3×10 −4 Pa and the deposition rate for Cr, SiO2 and Au was 0.5 Å/s.The main steps of the device fabrication are shown in Fig. S1.The magnetic anisotropies in the protected and gated regions were confirmed with polar MOKE measurements (Fig. S1).

MFM measurements
The MFM measurements were performed with a Bruker Dimension Icon Scanning Station mounted on a vibration-and sound-isolation table using tips coated with CoCr.To minimize the influence of the stray field from the MFM tip during the measurements, low-moment MFM tips were adopted.We repeated the MFM measurements and found that the magnetization in the nanomagnets remains unchanged, confirming that the MFM tips do not alter the magnetic configurations.For the voltage-control of the coupled nanomagnets, they were mount on a dedicated holder and connected to a source meter (Keithley 2400) with wire bonding.The MFM images were captured after employing the demagnetization protocol (Supplementary Information S2).All of the MFM measurements were performed at room temperature and under ambient conditions.

Electrical measurements
For electrical measurements, the magnetic films were patterned onto a 1.5 μm-wide Hall cross using electron-beam lithography and the coupled nanomagnet elements were located in the centre of the Hall cross.The devices were then connected to a source meter (Keithley 2400) and voltmeter (Keithley 2182) with wire bonding.All of the electrical measurements were performed at room temperature and under ambient conditions.

Micromagnetic simulations
To understand the mechanism of the AP/P coupling conversion, micromagnetic simulations were carried out with the MuMax3 code 42

S2. Demagnetization protocol
In order to obtain the low-energy magnetic configuration in an array of coupled nanomagnets, an oscillating magnetic field is applied perpendicular to the devices and the field amplitude is reduced over time (Fig. S2a).For this demagnetization protocol, the oscillating frequency of the magnetic field is 2 Hz (oscillation period t0 = 0.5 s) and its amplitude is linearly reduced from Hmax = 200 mT to zero with a constant step size of ∆HDemag = 0.167 mT.As shown in the MFM image of as-fabricated magnetic configuration of the Ising square-lattice element prior to applying magnetic fields (Fig. S2b), AFM-like domains are spontaneously formed, indicating the presence of AP coupling in the as-fabricated state.Following application of the demagnetization protocol, the size of the AFM-like domains increases, indicating that the array of coupled nanomagnets has been driven to a lower-energy state (Fig. S2c).Interestingly, on repeating the same demagnetization process, the AFM domain pattern changes, implying that the formation of AFM-like domains is stochastic.where ts, f0, Eb, k and T represent the time to switch, the attempt frequency (~10 9 Hz), the energy barrier to switching, the Boltzmann constant and temperature, respectively.The switching energy barrier can be determined using the Stoner-Wohlfarth model where: when the effective magnetic field Heff is parallel to the direction of the original magnetization and when the effective magnetic field Heff is antiparallel to the direction of the original magnetization.Here, the effective magnetic field Heff includes the external magnetic field Hext and the combined effect of the coupling with four nearest-neighbour sites: The magnetic moment m on the square-lattice site, the nearest-neighbour AP coupling strength J (-2.5 eV) and the anisotropy-induced switching energy barrier Esw (15.7 eV) are all taken from the experimental results.Here, J and Esw on each site is given by a Gaussian distribution to take into account the disorder in real devices.The magnetic configurations are obtained for different demagnetization step sizes ∆HDemag.The nearest-neighbour correlation <SiSi+1> is determined to evaluate how close the magnetic configuration is to the ground state.The AP ground state on the square lattice is well-defined, forming a "checkerboard" pattern with <SiSi+1> = -1.As shown in Fig. S3a, the magnetic configuration approaches the low-energy state on decreasing the demagnetization step size.
The demagnetization process can behave as a "thermal bath" that allows the array of coupled nanomagnets to relax into a low-energy configuration at an effective elevated temperature Teff 55,56 .We employ the Metropolis-Hastings algorithm to estimate Teff for our demagnetization protocol using the same coupling strength and distribution used in the macrospin model.The change of <SiSi+1> with respect to the effective temperature parameter of βJ (β = 1/kT) is shown in Fig. S3b.A transition in <SiSi+1> occurs around βJ ≈ 0.5, which agrees with the theoretical prediction of phase transition in the square-lattice Ising model: By comparing the values of <SiSi+1> obtained using the macrospin model and using the Metropolis-Hastings algorithm, we can estimate the effectiveness of demagnetization protocol to drive the coupled nanomagnets to their ground state (Fig. S3c).The decrease of the demagnetization step size effectively decreases the temperature of "thermal bath" that saturates at a certain temperature (ꞵJ ≈ 0.44) related to disorder in the coupled system.It also verifies the effectiveness of our demagnetization protocol with the experimental demagnetization step size of 0.167 mT being sufficient to realize the lowest-energy magnetic configuration.The protected regions are designed to have a high perpendicular magnetic anisotropy, which ensures that the magnetization is not perturbed by the stray field from the MFM magnetic tips during the measurements.This also means that the energy barrier for magnetization switching in the protected regions is higher than the thermal energy at room temperature.The coupled nanomagnet system is thus athermal and no thermally-active magnetization switching is observed during the experiments.Furthermore, the nearest-neighbour coupling strength is weaker than the energy barrier for switching the magnetization.Therefore, the voltagecontrolled change of the coupling strength cannot induce the spontaneous switching of the magnetization without applying the demagnetization protocol.
As shown in Fig. S4a, the array of coupled nanomagnets exhibits an AFM-like pattern following demagnetization.The device was then exposed to a negative voltage of -2.5 V for 120 min, converting the coupling from AP to P. The same magnetic configuration was observed in the subsequent MFM measurement indicating that the energy barrier for switching the magnetization is higher than the coupling strength and the thermal energy (Fig. S4b).In order to demonstrate that the nearest-neighbour coupling has switched from AP to P, the array was demagnetized again and an FM-like pattern was observed, confirming the change of coupling from AP to P (Fig. S4c).Therefore, experimentally, demagnetization of the array is essential to show the conversion of the voltage-controlled coupling.

S3. Details of macrospin and semi-micromagnetic model
In this section, we will first give a more detailed description of the macrospin model that was briefly introduced in the main text and shown in Fig. 3a.We then turn to a semimicromagnetic model to quantitatively estimate the coupling strength.
In the macrospin model, due to the strong OOP magnetic anisotropy, S1 and S2 can only  In Pt/Co, the antisymmetric exchange interaction is weaker than the symmetric exchange interaction, i.e. |Jex| > |Deff|.The energy curves for the AP and P configurations are shown in Fig. S5.By determining the difference in energy between AP and P alignment, we can obtain the strength of the coupling between S1 and S2.As discussed in the main text, if the gated region is strongly IP (Kg << 0), J ≈ Deff < 0, whereas if the gated region is strongly OOP (Kg >> 0), J ≈ Jex > 0, so providing an intuitive picture for the AP/P coupling conversion resulting from the Kg-mediated competition between symmetric and antisymmetric exchange interaction.
Despite the fact that it is possible to explain the AP/P coupling conversion with the macrospin model, the "effective" interaction terms of Jex and Deff in Eq. 3 are not clearly related to the material parameters in real devices.To get closer to a real physical system, a semimicromagnetic analytical model is developed on the basis of the macrospin model.In the semi-micromagnetic model, the total energy including exchange energy, anisotropy energy and DMI energy, can be written as: ( ) where A, K and D are the exchange energy constant, anisotropy constant and the DMI constant, respectively.A and D are constant throughout the magnetic regions, while K is different in the protected and gated regions.We denote the anisotropy constant within the gated region as Kg and as K0 for the protected region.
For simplicity, we assume the magnetization lies in xz plane and rotates linearly within domain walls.The structure of the basic element used for the model is shown schematically in Fig. S6a.
Substituting Eq.S9 into Eq.S8, we obtain the following expression for the energy: where S is the cross-sectional area.
When Kg < 0, and after the total energy minimization with respect to LDW, one finds: Since the gated region is narrow, we assume that LDW > wg and ( ) Taking LDW into account, the energy expression becomes ( ) When the OOP magnetic anisotropy in the gated region is relatively strong (Kg > 0), the domain wall will be fully located in the gated region i.e., LDW < wg.The expression of energy can then be written as: By minimizing the total energy with respect to LDW, we find: LDW into the energy density, we obtain: (ii) For P alignment (S1 = ↑ and S2 = ↑) (Fig. S6c), the boundary condition is According to the macrospin model, 0 2 π θ = ± when Kg << 0 and 0 0 θ = when Kg >> 0.
We first consider the case where Kg << 0 i.e., when the gated region is IP magnetized and 0 2 π θ = .Substituting Eq.S14 into Eq.S8, we obtain the following expression for the energy density: After minimizing the total energy with respect to LDW, we find: DW 0 2A L K π = .Again, since the gated region is narrow, we assume LDW > wg and ( ) LDW into the energy density expression, we obtain: ( ) We then consider the case of Kg >> 0 i.e., when the gated region is OOP magnetized and take 0 0 θ = .Substituting Eq.S14 into Eq.S8, we obtain the following expression for the energy: By determining the energy difference between the AP and P configurations, we obtain: and Comparing this with the results obtained from macrospin model, we obtain the relationship between physical parameters D and A, and the "effective" interaction terms Deff and Jex: The validity of these equations is confirmed by the good agreement of the coupling strength with that obtained from general micromagnetic simulations as described in main text.In particular, the magnitude of the AP coupling is given by: mT/2).In addition, we find that the P coupling strength has a g K dependence (Fig. 3b).Since the value of Kg tends to saturate at a certain value for VG < 0, the coupling strength J for the P coupling should have an upper limit.

S5. Interplay between DMI and magnetic anisotropy
In main text and Fig. 3, we present relationship between coupling strength and magnetic anisotropy in the gated region determined from micromagnetic simulations.The effective OOP magnetic anisotropy Keff used in the micromagnetic simulations is given by: where Ku and μ0 are the interfacial uniaxial magnetic anisotropy constant and the magnetic permeability of free space, respectively.
In Fig. S8, we present further results from the micromagnetic simulations, highlighting the interplay between DMI and magnetic anisotropy in the gated region.When Kg < 0, the energy of the system for AP and P alignment is almost the same in the absence of DMI, which supports the fact that DMI is responsible for AP coupling.With increasing DMI, the difference in energy for AP and P alignment increases when Kg < 0, indicating the enhancement of AP coupling.
This leads to an increase in the critical Kg where the coupling is converted from AP to P. When Kg > 0, the energy of the system for AP alignment surpasses that for P alignment, resulting in a P coupling.
In the main text, we only consider the VCMA effect for the voltage control of the coupled nanomagnets.However, it has been reported that the DMI strength can be modified with electric fields and that the DMI decreases with decreasing OOP magnetic anisotropy since these two effects share the similar origin of spin-orbit coupling [57][58][59][60] .This may be partially the reason for the slight difference between experimental (-2.5 eV) and calculated values (-3.0 eV) of the AP coupling strength.

S6. Reliability of <SiSi+1> obtained from different chips and positions
To illustrate the device-to-device reliability, the nearest-neighboring correlation function <SiSi+1> of the square lattice of 3 chips and 5 different devices per chip were measured as a function of the gate voltage (Fig. S9).The performance of the voltagecontrolled magnetic coupling on changing the gate voltage is found to be robust with the trend in <SiSi+1> reproduced within <0.1 standard deviation on the same chip.

S7. Effect of the dipolar interaction
Here, we determine the effect of the dipolar interaction in the Ising artificial spin ice.First, we estimate the contribution of the dipolar interaction for the basic element with two protected regions (Fig. 1a and 1b).A rough estimation of the energy of the dipolar coupling between the two protected regions can be obtained by considering two point-like dipoles placed at the centre of each element at a distance r from each other.In this case, the dipolar coupling Jdip is given by:  We now estimate the dipolar interaction in the extended Ising artificial square ices.For this, we consider the dipolar and exchange interactions in square lattice shown in Fig. S10.As the dipolar interaction is a long-range interaction, the energy associated with it needs to take into account the interactions from the surrounding macrospins.We evaluate the effect of the dipolar interactions by calculating the energy difference when flipping the central macrospin where ∆Edip and ∆Eex are the change in the dipolar and exchange energies on flipping the central spin.For simplicity, we consider the case where the surrounding macrospins are in the ground state.For AP coupling, the ground state has AFM order, and the change in energy due to the dipolar interaction on flipping the central macrospin is: where si and ri represent the orientation of the i th surrounding macrospins and the distance between the center and the i th macrospin, respectively.The dipolar energy is summed over all surrounding macrospins in a 100 × 100 square lattice similar to the experimental size.The dipolar interaction facilitates the formation of the AFM order.Due to the alternating up-down alignment of the magnetization for AP coupling, the energy change on flipping the central spin due to the dipolar interaction is small compared to ∆Eex ≈ 8J ≈ 19.6 eV.For P coupling, the ground state has FM order, and the energy change due to the dipolar interaction on flipping the central macrospin is: The dipolar interaction inhibits the formation of the FM order and, due to the uniform alignment of the macrospins in the P state, the energy difference induced by dipolar interaction becomes sizable.
In the experiment, we find that the strength of the AP and P coupling measured in the basic element is similar (Fig. 2d), while the correlation function <SiSi+1> for FM order is significantly smaller than that for AFM order in the extend structures (Fig. 4f).This could be due to the fact that the dipolar interaction becomes considerable in extended lattices and inhibits the formation of the FM order.To demonstrate a complex network with crossed connections, we present in Fig. S13a and S13b, an Ising-like nanomagnetic network including 5 vertices (Si; i=1…5), 6 magnetic connections (J12, J23, J34, J45, J35, J15; Jij <0) and 1 electrical connection (J14 <0).The nanomagnetic structure is similar to that shown in Fig. 6, which also contains two rings with an even and odd number of vertices.MTJs are fabricated on each spin vertex, and can be used to read and write the magnetization of the underlying Ising element (free layer) via the tunnel magnetoresistance effect and the spin transfer torque (STT) effect, respectively (Fig. S13b).

S8. Programmable coupling configuration in a four-spin chain
Each MTJ is addressed by a current Ii.When Ii is small, the STT effect is negligible and the magnetization of the Ising element (free layer) can be read via the tunnel magnetoresistance effect.When Ii is large, the STT effect tends to switch the magnetization of the Ising layer parallel or antiparallel to the reference layer, depending on the polarity of Ii.Assuming the magnetization in the reference layer to be↑, a positive (negative) Ii gives an effective magnetic field pointing ↑ (↓) whose strength is determined by the intensity of Ii.In addition, the ionic gate structures are fabricated on each connection region (shown in yellow in Fig. S13b) and the gate voltage Vij is used to tune the magnetic coupling Jij between the vertices Si and Sj.
Coupling J14 between vertices S1 and S4 is not possible in the 2D structure.Instead, we can couple vertex S1 and vertex S4 by applying electric currents through the MTJs on vertex S1 and vertex S4.The ground state of the Ising network is then obtained by applying the demagnetization protocol.The magnetization of S1 and S4 is read via the MTJ resistance with a small electric current.In addition, the electric currents are injected to couple S1 and S4 via the STT effect.Here, I14 is the intensity of the electric current corresponding to the coupling strength J14 and sign(Si) determines the polarity of the electric current.As J14 <0, I14 <0 and the electric currents I1 and I4 cause the vertices of S1 and S4 to be AP.For example, when S1 =↑, I4 =I14 <0 and S4 experiences an effective magnetic field pointing ↓, resulting in a current-induced AP coupling.This electrical process, including magnetization reading and electrical coupling, is continuously repeated throughout the demagnetization process to accomplish the integration of magnetic and electrical couplings.
In order to verify the effectiveness of the hybrid MTJ/Ising network structure, we performed a simulation with the same microspin model used for the square lattice.For the case of all magnetic couplings, the demagnetization process gives the four degenerate low energy spin states with approximately equal percentages as shown in Fig. S13c.In the presence of the electrical coupling J14, the electric currents I1 and I4 effectively couple the magnetization direction of S1 to that of S4.In the simulations, we set the time period of updating the electric currents to be 1 ms, which is faster than the response time of magnetic field in our experiment.
The STT-induced effective magnetic field is set to match the coupling strength of the magnetic coupling.After the demagnetization process, the simulation yields the doubly-degenerate low energy spin state with approximately equal percentages (Fig. S13d).
For a more general and complex Ising network, the vertex Si has Ni virtual couplings interacting with S i 1…S i Ni, and the electric current Ii is given by: Following the same principle, we can create a controlled-Majority gate, a functionally complete logic gate where any Boolean function can be implemented using a combination of Minority gates (Fig. S14d).The output of the controlled-Majority gate depends on the relative alignments of three inputs.As shown in Fig. S14e, if an electric voltage "1" ("0") is applied to the gate electrode, the coupling is set to be AP (P) and the direction of the output magnetization is opposite (equal) to that of the majority of the three input magnetizations, accomplishing the Minority (Majority) operation.Therefore, our logic scheme has the capability of dynamic reconfigurability, which increases the logic functionality of a device without increasing the number of logic gates and lead to a more compact logic chip.

Finding the solution to
the Max-Cut problem consists of maximizing the total weight of the edges between two mutually exclusive subsets of vertices.In our implementation, the nanomagnets represent the vertices of the graph, which are separated into two sets according to their magnetization, ↑ or ↓.The coupling strengths Jij corresponds to the weight of the edges wij in the graph.Solving the Max-Cut problem is equivalent to minimizing the energy of an Ising network with the same connections.Since we can electrically program the strength of each individual coupling, our nanomagnetic system can serve as a combinatorial optimization problem solver.In Figure6aand 6b, a schematic of an eleven-spin Ising network is shown, which represents a specific Max-Cut problem where each nanomagnet has either two or three connections.As shown in Fig.6c(upper panel), the demagnetized configuration of the Ising network whose couplings are all AP (Jij < 0) reveals the solution for the Max-Cut problem of a graph whose edge weights are all positive (wij > 0).In order to change the weight w34 to negative (w34 < 0), the corresponding coupling J34 is tuned to P (J34 > 0) by applying VG = -2.5 V for 90 min and the demagnetized magnetic configuration gives the new solution (Fig.6c, lower panel).

Figure 1 |
Figures and figure legends:

Figure 2 |
Figure 2 | Reversible conversion of AP/P magnetic coupling.a, Schematic of the Hall device used for electrical transport measurements.b, Magnetic hysteresis loops of a nanomagnet element on the Hall bar with AP (top) and P (bottom) coupling.The full hysteresis loops show three resistance levels corresponding to magnetic configurations of ↑↑, ↑↓ (↓↑) and ↓↓.The half hysteresis loops are measured after employing the demagnetization protocol and with the magnetic field starting from zero.c, Minor hysteresis loops with magnetic field starting from saturation.The magnitude of Hbias is indicated with the black dashed lines.The Hbias obtained from positive minor hysteresis loops (top 5 curves) is opposite to that obtained from negative minor hysteresis loops (bottom 5 curves).d, Evolution of Hbias obtained from positive and negative minor hysteresis loops with respect to the gate voltage.The error bars represent the uncertainty in the estimation of Hbias.Red and black dashed lines are guides to the eye.e, Percentage of AP magnetic configurations obtained after demagnetization (lower panel) with the gate voltage given in the upper panel.Each percentage is determined from 30 trials carried out on one device.Red-and blue-shaded regions highlight the application of positive and negative gate voltages, respectively.f, Stochastic behaviour of the demagnetized configuration for different coupling strengths.The labels I, II and III correspond to the states indicated in e.

Figure 3 |
Figure 3 | Macrospin model of AP/P coupling conversion and micromagnetic simulations of a nanomagnet element.a, Schematic of macrospin model consisting of three macrospins representing the magnetizations in two protected regions and one gated region.b, Total energy for the AP (blue) and P (red) alignment as a function of Kg determined from the micromagnetic simulations.The difference between the AP and P energies gives the coupling strength J (green).The AP coupling strength saturates at -3.0 eV, as indicated by the horizontal black dashed line, whereas the strength of P coupling increases with

Figure 4 |
Figure 4 | Voltage-controlled magnetic phase transition.a, Coloured SEM image (top) and corresponding schematic (bottom) of a 1D Ising-like chain structure.The widths of the protected and gated regions are 150 nm and 50 nm.b, MFM images of the chain element with AP (top) and P (bottom) coupling.c, Nearest-neighbour correlation function <SiSi+1> of AP and P coupling in a chain of 30 coupled regions as a function of gate width wg.The dimension of the gated region is indicated in the inset.The red and blue dashed lines are guides to the eye.d, SEM image of the 2D Ising-like square lattice structure.Part of the image is indicated in colour (left) with a corresponding schematic (right).e, MFM image sequence of the 2D Ising-like square lattice showing reversible magnetic phase transitions between AFM and FM order.f, <SiSi+1> as a function of gate voltage in the 2D square lattice.The states corresponding to the MFM image sequence (I to V) are indicated.The error bars represent the standard deviation of <SiSi+1> evaluated from five 15×15 lattices.The bright and dark areas in the nanomagnet regions in the MFM images correspond to ↑ and ↓ magnetization, respectively.Red-and blueshaded regions in the SEM images indicate the protected and gated regions.All the scale bars are 1 µm.

Figure 5 |
Figure 5 | Addressable control of magnetic coupling in a four-spin chain.a and b, Coloured SEM image and corresponding schematics of a programmable four-spin Ising chain.Red-and blue-shaded regions indicate the protected and gated regions, and yellow-shaded regions indicate the gate electrodes.c, Programming rules for the gate voltage and coupling.d and e, Coupling configurations for "111" (d) and "011" (e) programmed by applying the corresponding electric voltages.The blue and yellow connecting lines represent AP and P coupling, respectively.The percentages of AP alignment for the spin pairs S1|S2, S2|S3 and S3|S4 after demagnetization are shown (left), illustrating the programmed coupling configuration.Each percentage is obtained from the measurement of 32 elements.The MFM images of four element structures are shown with the bright and dark areas in the nanomagnet regions correspond to ↑ and ↓ magnetization, respectively, which is indicated with green and purple arrows (right).In order to guarantee the complete AP/P conversion, the gate voltages are applied for 90 min.All the scale bars are 500 nm.

Figure 6 |
Figure 6 | Programmable Ising network as hardware solver for Max-Cut problems.a, Schematic of a programmable Ising network based on coupled nanomagnets.Each coupling strength can be programmed by applying the corresponding electric voltage.b, Coloured SEM image and corresponding schematic of programmable eleven-vertex Ising network.In the SEM image, red-and blue-shaded regions indicate the protected and gated regions, while the yellowshaded region indicates the gate electrode.c, Solutions to Max-Cut problem obtained from MFM images of demagnetized devices for the cases when J34 is programmed to be AP (top) and P (bottom).The blue and yellow connecting lines in the schematics represent AP and P coupling.The black dashed line in each of the schematics indicates the cut lines separating vertices into two complementary sets (in green and purple), which is the solution to the Max-Cut problem with the corresponding weights.The bright and dark areas in the nanomagnet regions in the MFM images correspond to ↑ and ↓ magnetization, respectively, which is indicated with green and purple arrows.All the scale bars are 500 nm.

Figure S2 |
Figure S2 | Demagnetization protocol.a, Schematic showing the demagnetization procedure applied to the coupled nanomagnets.b, MFM image of the as-fabricated magnetic configuration of the Ising artificial spin ice prior to applying the demagnetization protocol.c, MFM images of the same area shown in b after applying the same demagnetization protocol.AFM-like domains are shaded in green and red.The bright and dark contrast in MFM images indicates nanomagnet regions with ↑ and ↓ magnetization, respectively.The scale bar is 1 µm.

Figure S3 |
Figure S3 | Simulation results of square-lattice Ising model.a, <SiSi+1> as a function of the demagnetization step size ∆HDemag.b, <SiSi+1> as a function of the effective temperature βJ obtained with the Metropolis-Hastings algorithm for the square-lattice Ising model.c, Effective temperature parameter βJ as a function of the demagnetization step size ∆HDemag.

Figure S4 |
Figure S4 | Demagnetized configurations following the electric gating.a, MFM image of the magnetic configuration of the as-fabricated square-lattice array following demagnetization.b, MFM image of the magnetic configuration in the same area after applying a negative voltage of -2.5 V for 120 min.c, MFM image of the magnetic configuration in the same area following a second demagnetization.AFM-like domains are shaded in green and purple in the zoomed-in regions of a and b, and FM-like domains are shaded in yellow and blue in the zoomed in region of c.The bright and dark areas in the nanomagnet regions in the MFM images correspond to ↑ and ↓ magnetization, respectively.The scale bars are 1 µm.

Figure S5 |
Figure S5 | Energies for AP and P alignment, and coupling strength as a function of Kg obtained from the macrospin model with Jex = 1.5 eV and Deff = -1 eV.

Figure S6 |
Figure S6 | Semi-micromagnetic model.a, Schematic of the basic element used for semimicromagnetic model.b and c, Schematics of magnetization tilt angle θ for AP (b) and P (c) alignment of the magnetization in the neighbouring protected regions.
Considering the continuous magnetization rotation and the symmetry of the structure, the magnetization at the centre can be either ← or →.However, due to the lefthanded chirality in Pt/Co, the magnetization at the centre prefers to be ← i.e., − .Hence the magnetization can be written as: ˆ[sin ,0,cos ] = .The magnetization can then be written as: ˆ[sin ,0,cos ]

Figure S7 |
Figure S7 | Snapshots of micromagnetic simulations for different Kg.The direction of the magnetization is indicated by the colour wheel, and white and black correspond to ↑ and ↓ magnetization, respectively.

Figure S8 |
Figure S8 | Energy of the system obtained from micromagnetic simulations for AP and P alignment as a function of Kg for different DMI values.

Figure S9 .
Figure S9.<SiSi+1> as a function of gate voltage in 15×15 square lattices obtained from 5 devices fabricated on 3 different chips.Red, green and blue colours indicate the results obtained from 3 different chips, while the different symbols indicate the results obtained from 5 different devices on the same chip.
with m = 3.0×10 -17 A•m 2 (for nanomagnet dimensions of 150 nm × 150 nm × 1.5 nm) and r = wp + wg = 200 nm, where wp and wg are the widths of the protected and gated region, respectively.This dipolar coupling is almost two orders of magnitude smaller than the measured coupling of 2.5 eV as well as the estimated exchange-induced coupling of 3.0 eV.Therefore, the effect of the dipolar interaction is negligible in the basic element.

Figure S10 |
Figure S10 | Schematics of magnetic configurations in the AP and P state used to estimate the effect of the dipolar coupling in an Ising artificial square ice.

Figure S11 |
Figure S11 | Programmable coupling configurations in a four-spin chain.a to h, All 2 3 = 8 coupling configurations that can be programmed using electric voltages.The applied voltages and corresponding coupling configurations, as well as one of the ground states, are shown.The blue and yellow connecting lines represent AP and P coupling, respectively.The probabilities of AP alignment for the pairs of S1|S2, S2|S3 and S3|S4 after demagnetization are shown (left), illustrating the programmed coupling configuration.Each percentage is obtained from the measurement of 32 elements.The MFM images of four selected element structures are shown with green and purple arrows indicating the magnetization of ↑ and ↓ respectively (right).In the MFM images, the bright and dark areas in the nanomagnet regions in MFM images correspond to ↑ and ↓ magnetization, respectively.In order to guarantee the AP/P conversion, the gate voltages are applied for 90 min.All the scale bars are 500 nm.

Figure S13 .
Figure S13.Hybrid MTJ/Ising network structure for solving a complex Max-Cut problem.a, Schematic of a 5-vertex Ising network.b, Schematic of the hybrid MTJ/Ising network structure.c and d, Simulation results of the ground state without (c) and with (d) electric coupling J14.The percentages of the ground spin state obtained from 100 simulation trials are indicated.

Figure S14 .
Figure S14.Reconfigurable nanomagnetic logic gates.a, Schematic of a controlled-NOT gate and corresponding logic circuit.b, MFM images of NOT operation for the AP coupling (top) and COPY operation for the P coupling (bottom).c, Truth table for the controlled-NOT gate.d, Schematic of a controlled-Majority gate and corresponding logic circuit.e, MFM images of Majority-NOT operation for the AP coupling (top) and Majority operation for the P coupling (bottom).In the schematics shown in a and d, red-and blue-shaded regions are the protected and gated regions, while yellow-shaded regions are the gate electrodes.The bright and dark areas in the nanomagnet regions in MFM images correspond to ↑ and ↓ magnetization, respectively.The blue and yellow lines in MFM images indicate the AP and P coupling.All scale bars are 500 nm.