Main

The dipolar interactions of a given spin with all of its nearest neighbours cannot be satisfied on a triangular lattice, resulting in a frustrated magnetic state with strong correlations and a local ordering principle, but no long-range order. Owing to its equivalence to the electrical charge distribution in water ice6, the materials are known as ‘spin ices’ and the local ordering scheme as the ‘ice rules’. Spin-ice materials such as Dy2Ti2O7 have been subject to an intense research effort7,8 and frustrated magnetism has evolved into a deeply interdisciplinary field, providing model systems for complex biological problems and a mathematical basis for the neural network algorithm from the Sherrington–Kirkpatrick model9.

A powerful way to understand spin ice is to consider the magnetic dipole as a positive and negative magnetic charge (±q) separated by one lattice spacing. The ice rule can then be described as the local minimization at each lattice site i of the total magnetic charge (Qi,qi). Predictions suggest that the magnetic properties can be fractionalized, with mobile excitations carrying magnetic charge, rather than spin, and their interactions being described by a magnetic Coulomb’s law2 (equation (1))

where V0 is the self-energy and ri j is the separation. Although these topological excitations are confined to the dipole lattice, and they are compatible with Maxwell’s equations10, their free magnetic charge character has led to the nomenclature magnetic monopole defects. Recent studies3,4,5,10 in rare-earth pyrochlores strongly suggest that monopole defects exist in bulk spin ice10.

Creating an odd number of intersecting dipoles, as in ‘kagome spin ice’11, is interesting because there is finite charge Q at each vertex, and a magnetic charge-ordered ground state of alternating Q=+q and Q=−q vertices (Fig. 1a and b) is predicted12,13. Figure 1c illustrates kagome ice-rule-defect (3-in or Q=+3q) creation. However a majority (diagonal), rather than a minority (vertical), bar flip would still create a magnetic charge carrier in the dipole strings, interchanging Q=+q and Q=−q vertices, but no ice-rule violation. In the bulk spin ices, the monopoles, magnetic charges and topological defects in the dipole string all refer to the same emergent excitations of the magnetic structure. In our artificial kagome ice, the magnetic charge and topological defects are conserved. The magnetic charge carrier can reside either at a vertex, in which case there is a complex magnetic charge distribution with a Q=±1 ice-rule state or a Q=±3 monopole ice-rule defect, or within a bar, in which case it is a transverse domain wall. Magnetization reversal does not require monopole formation.

Figure 1: Schematic of magnetic charge order and monopolar defect formation in artificial kagome spin ice.
figure 1

a, Schematic of a magnetic charge-ordered state of alternating Q=+q (red) and Q=−q (yellow) vertices of a 2D artificial spin-ice honeycomb nanoarray. b, We use an in-plane field (into the page) to split the manifold and select the saturated Ising state with parallel total vertex moments M=Σm. c, By flipping the central moment with an applied field out of the page, two ‘monopole’ defects are created.

Artificial spin-ice systems14,15 consist of regular arrays of magnetic squares15 or honeycomb lattices16, with dimensions in the 100–1,000 nm range, small enough to ensure single-domain behaviour but with sufficient magnetic volume for stable room-temperature ferromagnetism. These systems show ice-rule behaviour16 but direct observation of monopole defects is lacking. Here we have fabricated two-dimensional (2D) cobalt honeycomb (kagome spin ice) nanostructures and we create and image monopole defects on the spin-ice lattice, in situ in the magnetic force microscope (MFM).

Honeycombs were fabricated by electron-beam lithography and processing of sputtered 20 nm Co films (see the scanning electron micrograph in Fig. 2a). The bars are single domain, magnetized along the long (Ising) axis and their energy barrier to switching is large compared with both unpatterned Co and the thermal energy, kBT. The absence of thermalization makes the many degenerate ground-state configurations inaccessible. Previous experiments have obtained a low in-plane magnetization state by using a large out-of-plane conditioning field16 or by ‘shaking-out’ with a small oscillating in-plane field17,18. Shaking-out promotes dipole disorder, but does not always yield the ground-state even for the smallest sections of honeycomb19. Our protocol is, instead, to fully magnetize the array in one direction, giving a magnetized, magnetic charge-ordered state (as shown in Fig. 2b), and then sweep the magnetic field in the opposite direction (partial loop) to catch the magnetic structure in a mid-transition low-magnetization state, and bring the field to zero (partial reverse loop). Figure 2c shows the constructed magnetization field loop for our 2D array and field history paths (a–f). If we return directly to the same precise conditioning field, we observe no defect movement; however, if we cycle the full hysteresis loop and then repeat the conditioning protocol we do observe differences in the images, consistent with the lattice being frozen in the absence of a field greater than the last applied switching field and a degree of randomness in the switching process.

Figure 2: Characterization of a honeycomb array.
figure 2

a, Scanning electron micrograph of cobalt honeycomb on Si. Bar dimensions a=1,000 nm, width=100 nm, thickness=20 nm. Shape anisotropy causes the bar to be single domain, with moment, m, along the long (Ising) axis. Each Ising spin acts as a bar magnet with magnetic charge ±qm=m/a at each end. b, MFM image in zero field after saturation in-plane (at H=0 on the black curve in c). Inset: Enlarged image showing the Q=−1 plaquette and all unscreened neighbours (Q=±1). c, Field in-plane normalized MH loop of the honeycomb array. The sweep direction is negative (black) and positive (red). The MH loop of the array has been calculated from anisotropic magnetoresistance data (Supplementary Fig. S1). Points a–f indicate the field conditioning profiles undertaken before collecting the MFM images in Fig. 4a–f. The black line is the simulated hysteresis loop taking a normal distribution of bars with mean HC=60 mT, σ=10 mT.

In Fig. 3a we show the MFM image, and in Fig. 3b a cross-section when a monopole has been created. Here the (positive) phase shift is a factor of three greater than at one of the ice-rule vertices, interpreted as a Q=+3q (ice-rule forbidden) vertex. In continuous honeycombs that obey the ice rules it is impossible to uniquely assign the dipole moments from MFM (ref. 17) as the image is not sensitive to flips of infinite chains and closed loops of head–tail dipoles. However, a Q=±3q state is unique and allows exact solution locally.

Figure 3: The creation of monopole defects.
figure 3

a,b, MFM image (a) and linescan (b) in zero field. The change in phase of the tip, proportional to d2Bz/dz2, is approximately three times greater at the ‘bright’ yellow vertex, corresponding to a magnetic charge Q=+3q.

Incremental (0.6 mT in Fig. 4a–f) steps in the conditioning field profile drive the defects from trap to trap, as shown in our sequential images and schematic pathways (each colour indicates a separate magnetic charge pathway). The kagome geometry prevents the pathways crossing in a single switching transition. The conserved magnetic charge carriers have ΔQ=±2q and transverse domain walls are expected in our cobalt nanowires20. The two possible ΔQ=+2q defects (Q=+q,|M|=2 m and Q=+3q,M=0) are topologically equivalent and similarly the two ΔQ=−2q defects, whereas at any given site we image a total magnetic charge Q=(±q±2q). On reversal of the magnetization in the array structure, the magnetic charge carrier is conserved and its motion down the dipole string is field directed and trackable, but monopole defects appear only when the magnetic charge carrier is trapped at a vertex site with like charge despite strong magnetic Coulomb’s law repulsion. This is a rare random event and therefore one described by Poisson statistics.

Figure 4: The movement of magnetic charge.
figure 4

MFM images in zero field after the labelled conditioning fields. a, Two 3-in monopole defects with Q=+3q form (bright yellow spots, yellow dots), with strings of head–tail spins (blue and green) to opposite magnetic charges ΔQ=−2q (red dots). Further bar flips are required to make the schematic (right) match the observed data; the grey arrows indicate a (non-unique) trial solution. b, The ΔQ=+2q magnetic charges hop to the right, changing from Q=+3q to Q=+q. cf, Another ΔQ=−2q magnetic charge appears, tracing its own (orange) string to the left, until it is blocked by the blue string of the ΔQ=+2q carrier.

To investigate this in detail, we create a randomly disordered model system. We treat each bar i as an Ising spin and assign a pseudorandom coercive field HC i, within a normal distribution with expectation value, 〈HC〉=60 mT, and a standard deviation, σ=10 mT, taken from the experimental hysteresis loop in Fig. 2c. Averaging over 1,000 such pseudorandom configurations, we find that purely random disorder does not capture the fine details of the hysteresis curve (Fig. 2c) and bears no resemblance to the MFM images (Supplementary Movie). By making the defects both unfavourable to form and unstable with a tunable ice-rule modification to the coercive field (HI) with Hflip,i=(HCi+CHIDHI), where C and D are the number of monopoles created and destroyed respectively on switching the bar, one can tune the monopole concentration by selecting values of HI and σ. In fact, we observe that the ratio HI/σ is the sole controlling parameter in this respect. The simulated hysteresis loops are compared with the data in Fig. 5a and a significantly better fit is obtained with HI/σ1. The predicted mean number of monopoles (N(obs)) at each field is plotted in Fig. 5b. For HI=0 the mean is simply the configurational probability (P=1/8) multiplied by the number of vertices (507). A statistical analysis of the individual frames of the data in the Supplementary Movie indicates that the mean N(obs) for 500 vertices is somewhere between two and four; therefore, P0.01. Figure 5d shows the experimental N(obs) in each movie frame with the simulated mean N(obs) curves with HI=0.75σ, 0.9σ and σ. With these values the model also reproduces the observed hopping characteristics of the magnetic charge carriers.

Figure 5: Simulation of data with a model of the normal distribution of bar coercive fields (HC: mean HC=60 mT,σ=10 mT) plus a tunable ice-rule modification to the coercive field (HI).
figure 5

If switching a given bar would create a Q=±3q state then Hflip=(HC+CHIDHI). Field-driven switching was simulated for 1,000 pseudorandom coercive field distributions on a lattice of 507 vertices (to match the frame size in the Supplementary Movie) with different values of HI. a, Section of hysteresis loop data (symbols) and mean simulated hysteresis loops for selected HI (solid lines) and for D=0,HI=2σ (dashed line). b, Mean number of simulated Q=+3 defects (N(obs)) versus field; there is no asymmetry in the populations of Q=+3 and Q=−3 defects. For HI=0 the mean is P=1/8 multiplied by the number of vertices (507). c, Histogram of the number of times a given N(obs) was found for both Q=+3 and Q=−3 defects; all frames in the Supplementary Movie were counted. Inset: (black, left axis) Magnetic Coulomb potential of a domain wall moving through a single vertex with QDW=2q,Qv,max=1.3q at each neck22 (0.05 μm from vertex centre) and a domain-wall width of 0.1 μm. (Blue, right axis) Schematic of the ice-rule energy EI profile in the simulations. (Olive green dashed line, right axis) Schematic of the ice-rule energy EI profile in D=0 simulations. d, Measured N(obs) versus field (±5% inhomogeneity) for each frame from the Supplementary Movie (columns) and simulated mean N(obs) curves with HI=0.75σ , 0.9σ and σ from b replotted.

We now consider the pinning required to form energetically unfavourable monopoles. Domain-wall spectroscopy in permalloy nanowires demonstrated that remote magnetic charge provides effective domain-wall pinning21 only out to 100 nm and revealed a local pinning potential well within the magnetic Coulomb barrier at an isolated permalloy Q=+1 T-vertex22. We conclude that local vertex effects dominate and the inset of Fig. 5c shows the schematic of the 1D local magnetic Coulomb’s pinning potential8 and the approximation used in our simulations.

Finally, the creation, movement, local site stabilization and detection of such quasi-particles in 2D systems may lead to alternative approaches in data storage and magnetic logic, potentially even to Sherrington–Kirkpatrick-type neural-net hardware, particularly given the electrical connectivity of the honeycomb and the research effort in devices based on domain-wall motion.

Methods

Fabrication.

Samples were fabricated by sputtering 20-nm-thick Co thin films followed by standard electron-beam lithography and processing. The lateral extent of the honeycomb is 100 μm×100 μm with cobalt electrodes processed in the same write step. The growth rate was calibrated and the thickness confirmed by vibrating sample magnetometry of the unprocessed film.

Magnetic force microscopy.

MFM was carried out with a Digital Instruments 3100 series model and standard Veeco MESP tips were used at room temperature. The system has been fitted with a custom-made electromagnet. MFM was carried out at remanence after the application of successive in-plane magnetic fields. Previous experiments have obtained a low-magnetization state on the honeycomb lattice by ‘shaking out’ with a small oscillating field17,18. Our protocol is to fully magnetize the array in one direction and then sweep the magnetic field in the opposite direction to catch the magnetic structure in a mid-transition low-magnetization state. During the application of the magnetic field, the tip was withdrawn above the pole pieces of the electromagnet to ensure the tip magnetization was not perturbed by the external field. At this height, the stray field from the electromagnet was <0.005 T, which is an order of magnitude lower than the coercivity of the tip. MFM images of 25 μm×25 μm in size were taken of a central area of the sample. Images were taken on the same area many times, to ensure that the stray field from the tip was not switching the nanostructure being studied. In the MFM technique the change in phase of the tip, with respect to the cantilever oscillation, is proportional to d2Bz/dz2. Repulsive (attractive) interactions lead to a negative (positive) phase shift; in our images this corresponds to dark (light) and a negative (positive) net magnetic charge at the site.

Modelling.

Simulations were carried out on a kagome lattice of Ising spins of dimensions 30 bars by 36 bars (similar to the experimental data). Vertices at the edges were not counted, giving 507 active vertices. The coercive field of each bar was set using a Gaussian-distributed pseudorandom number Hc. The effective flipping field of each bar Hflip is thus determined by Hc and two further terms reflecting the energy cost of creating and destroying a Q=±3q state. Thus, Hflip=(Hc+CHcreateDHdestroy), where C and D are the number of Q=±3q states created and destroyed respectively on switching the bar (either 0 or 1). The array is initially negatively polarized and the values of Hflip are calculated for each bar. The bar with the lowest value is allowed to flip to positive polarity and then the lowest remaining value of Hflip is determined as before. In each case the applied magnetic field is recorded as the value of Hflip unless this is less than the previous value, in which case the magnetic field remains the same. As such, multiple flips at single values of applied field (cascades) are possible. This process is repeated until all of the bars have been flipped and then reversed to generate the downward cycle. The entire procedure was repeated 1,000 times and the results averaged. We find that the mean populations and lifetimes of Q=+3 and Q=−3 defects are symmetric and the observed Q=+3 population can be simulated with just a penalty to create an ice-rule defect (finite Hcreate), although these defects are generally too long-lived. The mean lifetime is shortened by introducing finite Hdestroy and a good agreement with the observed data is found when Hdestroy=Hcreate=HI; this is the situation described in the text. Independently, the Hcreate and Hdestroy terms have an almost identical effect on the mean population, and both allow the occurrence of cascades (which are completely absent if Hcreate=Hdestroy=0).