Introduction

The concept of magnetic monopole transport within a condensed matter setting has captivated scientists, allowing established theory1 to become an experimental realization2,3,4 within the bulk pyrochlore systems known as spin-ice5. In these three-dimensional (3D) systems, rare earth spins are located upon corner-sharing tetrahedra, and energy minimisation yields a local ordering principle known as the ice-rule, where two spins point into the centre of a tetrahedron and two spins point out. Representing each spin as a dimer, consisting of two equal and opposite magnetic charges (\(\pm q\)), is a powerful means to understand the physics of spin-ice5. Using this description, known as the dumbbell model1, the ice-rule is a result of charge minimisation, yielding a net magnetic charge of zero in the tetrahedra centre (\(Q={\sum }_{i}{q}_{i}=0\)). Then the simplest excitation within the manifold produces a pair of magnetic charges (\({\sum }_{i}{q}_{i}=\pm 2q\)) which, once created, can propagate thermally and only at an energy cost equivalent to a magnetic analogue of Coulomb’s law. The energy scale for the production of monopoles upon the spin-ice lattice is controlled by the chemical potential (μ), which is governed by properties intrinsic to the material such as lattice constant and magnetic moment6. Canonical spin-ice materials have a chemical potential that places them in a weakly correlated regime where only a small fraction of bound monopole–antimonopole pairs are found. Recent theoretical work has studied the ordering of magnetic charges upon cleaved spin-ice surfaces, perpendicular to the [001] direction7. In such systems, the orphan bonds upon the surface are found to order in either a magnetic charge crystal or magnetic charge vacuum, depending upon the scales of exchange and dipolar energies7. Experimental studies are now hinting at the presence of a surface-driven phase transition8 but the transport of magnetic charge across such surfaces has not been considered previously.

The arrangement of magnetic nanowires into two-dimensional lattices has recently shown to be a powerful means to explore the physics of frustration and associated emergent physics. These artificial spin-ice (ASI) systems9,10,11,12,13,14,15, where each magnetic nanowire behaves as an effective Ising spin, have recently yielded an experimental realisation of the square ice model16 and have also been used to study the thermal dynamics of monopoles in the context of Debye–Hückel theory17. Controlled formation of magnetic charge is an exotic means to realise advanced multistate memory devices. Such concepts have been shown in simple 2D lattices using magnetic force microscopy (MFM)18. The extension of artificial spin-ice into true 3D lattices that capture the exact underlying geometry of bulk systems is paradigm-shifting, allowing the exploration of ground state ordering and magnetic charge formation in the bulk as well as upon the surface. The production of 3DASI systems harbouring magnetic charge also allows marriage with advanced racetrack device concepts19,20.

In this study, we use state-of-the-art 3D nanofabrication and processing in order to realise a 3DASI in a diamond-bond 3D lattice geometry, producing an artificial experimental analogue of the originally conceived dumbbell model1. MFM is then harnessed to image the formation and propagation of magnetic charge upon the 3D nanowire lattice.

Results

Figure 1a shows a schematic of the 3DASI, which is composed of four distinct layers, labelled by colour. The system is fabricated by using two-photon lithography21,22,23,24 to define a polymer lattice in a diamond-bond geometry, upon which 50 nm Ni81Fe19 is evaporated (see ‘Methods’ for further details). This yields NiFe nanowires within a diamond lattice geometry as shown previously24. Each nanowire has a crescent-shaped cross-section (Fig. 1a inset), is single domain and exhibits Ising-like behaviour24. The L1 layer which is coloured red, is the upper surface termination and consists of an alternating sequence of coordination two (bipods) and coordination four vertices (tetrapods). The L2 and L3 layers, coloured blue and green, respectively, are ice-like with only vertices of coordination four. Finally, the L4 layer, coloured grey, is the lower surface termination of the lattice and consists of vertices which alternate between coordination two and coordination four.

Fig. 1: A 3D artificial spin-ice.
figure 1

a Schematic of a 3D artificial spin-ice system. The surface L1 layer (red) consists of an alternating sequence of coordination two and coordination four vertices. Below this, the L2 (blue) and L3 (green) layers can be seen. Within these layers, only vertices of coordination four are present. The L4 layer (grey) is the lower surface termination which again has an alternating sequence of coordination two and coordination four vertices. C2 and C4 label coordination two and coordination four vertices, respectively. Inset: Cross-section of Ni81Fe19 (grey) upon the polymer scaffold (yellow). b A false colour scanning electron microscopy image of the 3D artificial spin-ice lattice. Scale bar is 20 µm. c A false colour scanning electron microscopy image showing the L1 (red) and L2 (blue) sub-lattices, viewed at a 45° tilt with respect to the substrate plane. Scale bar is 1 µm. d Atomic force microscopy image of the 3D artificial spin-ice system. Scale bar is 2 µm. Coordinate system for field application is shown in top-right image. e Possibilities for creating magnetic charge upon L1. f The possible states and associated magnetic charge that can be realised at vertices of coordination two and coordination four.

Full size image

The overall array size is approximately 50 μm × 50 μm × 10 μm as seen in the scanning electron microscopy (SEM) image (Fig. 1b). Analysis of SEM data (see Supplementary Fig. 1) indicates the long axis of L1 wires is orientated at θ = (33.11 ± 2.94)° from the substrate plane, matching within error the angle of 35.25° which is expected for an idealised diamond-bond geometry5. A higher magnification image, clearly showing the L1 (red) and L2 layers (blue) can be found in Fig. 1c. The topography of the upper three layers can be measured using atomic force microscopy (AFM) as shown in Fig. 1d. The coordinate system used to define field directions is also shown in Fig. 1d.

The surface of the 3DASI lattice, which is composed of coordination two vertices upon L1 and coordination four vertices upon L2, provides interesting possibilities with respect to magnetic charge transport. In Fig. 1e, we illustrate how magnetic charge propagates along the L1 layer. Starting with a saturated state, applying a magnetic field above a critical value along the unit vector (1,−1,0) leads to the nucleation of a domain wall (DW) (Fig. 1e, top-left) which carries a mobile magnetic charge of magnitude ±2q. When reaching the L1–L2 junction (Fig. 1e, top-right), the effective vertex magnetic charge becomes Q= + 2q. A further increment in magnetic field leads to the L1–L2 junction emitting another DW (Fig. 1e, bottom-left) and when this wire is fully switched a surface magnetic charge state of Q= + 2q is realized (Fig. 1e, bottom-right). Note that a field applied in either direction along [\(110\)] with a projection along the L2 sub-lattice produces only magnetic charges at four-way junctions (see Supplementary Fig. 2). Overall, the 3DASI surface can realise effective magnetic charge magnitudes of ±4q, ±2q and 0 as summarised in Fig. 1f.

Imaging the magnetic configuration of a 3DASI

MFM is a convenient method to deduce the magnetization configuration of the 3DASI during field-driven experiments. This imaging technique is sensitive to the second derivative of the stray field with respect to z (d2Hz/dz2) which makes it ideal for imaging magnetic charge25 upon the 3DASI lattice. In the present study, we focus upon the field-driven transport of magnetic charge upon the L1 and L2 layers. The volume of the individual nanowires is sufficiently high that the 3DASI system is frozen at room temperature and thus thermal energies are negligible when compared to the energy required to switch a wire.

It is initially insightful to first study the simplest scenarios where each sub-lattice is saturated. Optical magnetometry (see Supplementary Fig. 3) indicates 30 mT is well above the saturating field for each sub-lattice. Figure 2a presents an MFM image, taken at remanence following a H = 30 mT in-plane magnetic field, first applied along unit vector (1,−1,0) and subsequently along unit vector (−1,−1,0). Masks are placed over void regions to guide the eye to signal originating from L1 and L2. Unmasked data is provided in Supplementary Fig. 4. Every L1–L2 vertex within the array is seen to have identical contrast. A magnified example of the contrast associated with an individual L1–L2 vertex is also shown in Fig. 2b, top. Here, the dashed lines separate individual islands and a corresponding 3D pictorial representation of the magnetic configuration is shown in Fig. 2b, bottom. With our choice of tip magnetisation, the bright yellow lobes indicate a positive phase associated with the stray field at magnetisation tail while bright red lobes indicate a negative phase associated with the stray field at magnetisation head. Focusing first upon the L1 nanowires, one can see lobes of strong positive contrast at the upper left of the nanowires and negative contrast in the lower right of the nanowires. Now focusing upon L2, strong positive contrast is seen in top right of nanowires, with negative contrast seen in bottom left. Overall, the vertex configuration is consistent with a type 2 ice-rule configuration produced by the applied field protocol. We note that near the bottom left of the L2 nanowires, faint positive contrast is seen (labelled A). A previous investigation, which took images in reversed tip configurations identified this as an artifact24, due to the abrupt upwards change in topography experienced by the tip at this point. Since the signal originating from the artifact is approximately a factor of two smaller than the signal originating from magnetic contrast, its presence does not impede analysis of the magnetic configuration.

Fig. 2: Imaging the saturated states in a 3DASI.
figure 2

a An MFM image taken at remanence after application of saturating fields along the unit vectors (1,−1,0) and (−1,−1,0). A coordination two, surface vertex is highlighted in pink, while a coordination four vertex at the intersection of L1 and L2 is highlighted in red. The scale bar represents 2 µm. b Top: Magnified example of the MFM contrast seen associated with L1–L2 junctions as seen in (a). Arrows are coloured by the local in-plane magnetization components. Dashed lines represent the separation between individual islands, with black for upper L1 wires and blue for lower L2 wires. Bottom: 3D Pictorial representation of the magnetization and MFM contrast. c MFM image taken at remanence after a further saturating field is now applied along unit vector (−1,1,0). d Top: Magnified example of the MFM contrast seen associated with L1–L2 junctions as seen in (c). Bottom: 3D Pictorial representation of the magnetization and MFM contrast. e MFM image taken at remanence after a further saturating field is now applied along unit vectors (1,−1,0) and (1,1,0). f Top: Magnified example of the MFM contrast seen associated with L1–L2 junctions as seen in (e). Bottom: 3D Pictorial representation of the magnetization and MFM contrast.

Full size image

To demonstrate that each sub-lattice can reverse independently, we now take images after saturating fields along different principal axes. Figure 2c shows the large scale MFM image taken at remanence after a saturating field along unit vector (−1,1,0). It is clear that contrast upon L1 wires have inverted. Further inspection of the magnified example (Fig. 2d) clearly shows the lobes of contrast upon L1 have indeed inverted showing the magnetization here has switched. This is also demonstrated pictorially in Fig. 2d, bottom. The contrast upon L2 is found to be unchanged, as expected. The system was then returned to the initial state (Fig. 2a) before a saturating field was applied along the unit vector (1,1,0). Examination of Fig. 2e now shows contrast upon every L2 nanowire has changed. Close inspection of Fig. 2f now shows stronger positive contrast in bottom left and strong negative contrast in top right, suggesting the wires have switched, as shown pictorially in Fig. 2f, bottom. Overall, these results provide confirmation that L1–L2 vertices corresponding to saturated states can be identified. Our previous work24 suggests that faint contrast is also expected at the top of L1 coordination two vertices (black dashed line in Fig. 2b, d, f) and at mid points upon L2, close to the L2–L3 junction (blue dashed line in Fig. 2b, d, f). Such contrast is expected even for uniformly magnetized states, due to a change in sign of Mz at the vertex. Upon L1, the effect of this is to smear out the edge contrast, such that fainter contrast of lower magnitude is seen at the L1 coordination two vertex. At the L2–L3 vertex, faint contrast is also seen, but we note that this is not currently sufficient to determine the magnetic state of the L3 layer.

With this fundamental understanding we next sought to understand the magnetic configuration of vertex states observed during the switching process to determine if monopole-excitations can be identified and tracked. Figure 3a shows an MFM image following a saturating field along the unit vector (1,−1,0) and subsequent 9.5 mT field along the unit vector (−1,1,0). Optical magnetometry indicates this is within the field range that switching is expected upon L1. A vector map of the magnetic configuration (Fig. 3b) has been produced through observations of the MFM contrast associated with each L1–L2 vertex as well as the surrounding wires. Unmasked images can be found in Supplementary Fig. 5. We note that there are multiple independent means to confirm the presence of a monopole. Firstly, contrast near the L1–L2 vertex is an excellent indication. If three of the four wires have contrast of the same sign, this is a monopole state. This can be further confirmed by then checking contrast upon the opposite ends of the wires. Finally, since the magnetic charge upon the wire ends closest to L1–L2 smears over the vertex area, the absolute magnitude of the phase is increased when compared to an ice-rule state. We have used all three criteria simultaneously to identify monopoles at the L1–L2 vertex. Interestingly, so long as a well-defined field protocol is used, it is also possible to infer the presence of monopoles at the L2–L3 vertex. Here, so long as L3 has been saturated, we expect this sub-lattice to be uniformly magnetised. However, if the extremities of two adjacent L2 nanowires both have positive or negative contrast, a monopole is implied at the L2–L3 vertex.

Fig. 3: Identification of monopole-excitations.
figure 3

a MFM image taken at remanence following a saturating field along the unit vector (1,−1,0) and subsequent 9.5 mT field along the unit vector (−1,1,0). b Associated vector map illustrating the magnetic configuration, monopole-excitations are annotated with bright yellow (Q = −2q) and red (Q = +2q) lobes. Each island represents a bipod, coloured with the local in-plane magnetization, as determined by key. c MFM image taken at remanence following a saturating field along unit vector (−1,−1,0) and subsequent 8.0 mT field applied along (1,1,0). The scale bars upon a, c represents 2 µm. d Associated vector map illustrating the magnetic configuration, and presence of monopole-excitations. e–h Magnified examples of the MFM contrast associated with L1–L2 junctions where Q = ±2q.

Full size image

In Fig. 3a, b, every L1–L2 vertex in the observed area resembles one of the patterns seen in Fig. 2b and d, with two exceptions. These are two monopole-excitations, each with a charge of Q = −2q, readily identified due to the enhanced MFM signal, which is a factor of 2 greater than the corresponding ice-rule state. Furthermore, the signal associated with the L1 wires on either side of the monopoles is clearly seen to oppose, whereas the L2 wires are identical and so must be aligned. Figure 3c, d shows a similar intermediate state following a saturating field applied along unit vector (−1,−1,0) and subsequent 8.0-mT field applied along (1,1,0). This allows intermediate states to be probed upon the L2 layer. Here, 9 monopoles are identified through observations of the contrast associated with each L1–L2 junction, as well as the surrounding wires. Figure 3e–h shows magnified examples of monopole-excitations with Q = ±2q, in each case, one pair of colinear wires exhibits opposing contrast with respect to one another, while the other pair of colinear wires show matching patterns of contrast. We note that for both intermediate states (Fig. 3a, c), the sub-lattice that extends along the field direction is effectively demagnetized (M < 0.1MS), so it is intriguing that a vast difference in the density of monopole-excitations is seen between the two images.

Tracking monopole propagation on the 3D lattice

To form a more complete understanding of the monopole behaviour on the surface, we now measure the detailed switching between two saturated states, taking images at 0.25-mT intervals. To do this we carry out direct observations of the reversal sequences for the L1 and L2 sub-lattices. Figure 4 shows vector maps representing snapshots (full MFM data can be found in Supplementary Figs. 6–9) of the switching process for the upper two layers of the lattice. Here each island corresponds to a bipod on the lattice, as defined in Fig. 1d. Each image contains ~70 wires on L1 and 70 wires on L2, only counting those where the majority of the wire is within the measured area. Analysis herein considers wires within this \(8\times 8\,\upmu {\mathrm{{m}}}^{2}\) measured region, this is due to a compromise between size of the observed area and data acquisition time. Figure 4a illustrates the array after application of 8.5 mT along the unit vector (−1,1,0). This field magnitude yields the first evidence of switching along this direction. Though much of the array remains saturated, six wires (three bipods) have switched yielding two monopole states, each with charge −2q (Monopoles 1, 2). In both cases, the monopoles are found at the intersection between L1 and L2. Further field increments yield additional chains of wires switching (Fig. 4b–d), with a further two negative monopoles (monopoles 3, 4) residing at the L1–L2 junction, after which L1 reaches saturation within the sampled area (Fig. 4e).

Fig. 4: Direct imaging of magnetic charge upon a 3D artificial spin-ice system.
figure 4

a–e Vector maps illustrating the magnetisation configuration and associated monopole-excitations in five snapshots during a reversal sequence upon the L1 sub-lattice. Here a saturating field was first applied along the unit vector (1,−1,0) after which a field of 8 mT was applied along the unit vector (−1,1,0). Successive images were then captured at remanence following 0.25-mT increments. Each island represents a bipod, coloured with the local in-plane magnetization, as determined by key. Each monopole excitation is assigned a unique index to track propagation between images. f–j Vector maps illustrating an equivalent reversal of the L2 sub-lattice. Here the samples were first saturated along the unit vector (−1,−1,0) after which a field of 6.50 mT was applied along (1,1,0). Successive images were then captured at 0.25-mT increments. Full datasets, including raw MFM images can be found in Supplementary Figs. 6 and 7.

Full size image

Figure 4f illustrates the measured region after the array had been saturated along the unit vector (−1,−1,0) and a field of 6.75 mT applied in (1,1,0). Eight monopoles can be immediately seen (monopoles 5–12), all of which seem to have appeared in pairs of ±2q. Here, five monopoles reside upon the L1–L2 junctions, while the remaining two reside upon L2–L3 junctions. Additional field increments lead to the creation of further monopoles (monopoles 12–18), while others move along the L2 nanowires or propagate out of the measured area (Fig. 4g–j).

The differences in monopole formation upon the L1 and L2 sub-lattices is striking. Application of an external field with component along L1 yields few uncorrelated magnetic charges (see Supplementary Fig. 10a) within the measured region, which seem to only be observed within a narrow field window (8–10.5 mT). We note that while this yields a net charge locally in the measured area, charge neutrality is expected across the full lattice. Analysis of the switching also shows a distinct absence of magnetic charges upon surface vertices with coordination two. On the contrary, the L2 switching leads to nucleation of many correlated pairs yielding almost equal numbers of positive and negative magnetic charges (see Supplementary Fig. 10b), meaning the net charge within the measured area is close to zero throughout the field range (see Supplementary Fig. 10c). The magnetic charges are also formed at a lower field (6.5 mT) for the L2 sub-lattice and remain for a wider field range (6.5–10.75 mT).

Modelling the 3DASI system

Calculating the total energy density of every possible vertex state, within a micromagnetic framework (see Supplementary Fig. 11) is an insightful exercise and provides some initial understanding of the system. Here it can be seen that the energy density to create a magnetic charge upon a coordination two, surface vertex is 3.2 times higher than that of a monopole at a coordination four vertex suggesting surface charges will be very unfavourable. To understand the significance of this within the context of switching the entire array, we carry out Monte–Carlo (MC) simulations based upon a compass needle model (see ‘Methods’). This is carried out for varying surface energetics factor (α) and quenched disorder arising from fabrication (\({d}_{i}\), see ‘Methods’). A disorder of \({d}_{i}\)  = 30% showed good agreement with switching field distributions in experimental data. The surface energetics factor (α) scales the energy required to produce a monopole upon the coordination two vertex, when compared to a coordination four vertex. A series of simulations with varying α are shown in Supplementary Fig. 12. Simulations which considered degenerate monopole surface energetics (α = 1, Supplementary Fig. 12) with a field applied along projection of L1 (−1,1,0) showed the presence of magnetic charges upon surface coordination two vertices and also short Dirac strings, in contrast to experimental data. Increasing the surface energetics factor to the value calculated in finite element simulations (α = 3.2), now reduces the number of magnetic charges seen upon surface coordination two vertices but Dirac string lengths are still shorter than seen in experimental results.

Figure 5a–c shows the results of MC simulations performed with enhanced surface energetics (α = 6.4) for field applied along the unit vector (−1,1,0). Upon the threshold of switching (Fig. 5b), chains of islands switch upon the L1 sub-lattice producing uncorrelated monopoles and long Dirac strings as seen in the experimental data before the majority of the array becomes saturated (Fig. 5c). Critically, charges upon surface coordination two vertices are now very rare, which is in agreement with experiment. Figure 5d–f shows MC simulations for the field aligned along unit vector (1,1,0). Here, a low field immediately produces large numbers of correlated monopole–antimonopole pairs (Fig. 5e), separated by a single lattice spacing, closely aligned with the experimental data. Figure 5g summarises the simulation results by showing the fraction of excited states obtained upon L1 and L2, showing excellent qualitative agreement with the experimental data, presented in Supplementary Fig. 10.

Fig. 5: Simulating the monopole dynamics upon a 3D artificial spin-ice.
figure 5

a–c Arrow maps showing magnetisation configuration with field applied along unit vector (−1,1,0). Arrows on L1 and L2 have black borders, arrows on L3 and L4 are borderless. d–f Arrow maps for field applied along (1,1,0) resulting in highly correlated monopole pairs. g Fraction of excited states during the L1 (blue) and L2 (red) reversal sequences.

Full size image

Discussion

As in all ferromagnetic materials, the 3DASI studied here passes through a field-driven state whereby the component along the field is effectively demagnetized. It is interesting to identify two main ways that this can be achieved in this 3D nanostructured system. The first possibility is that of local demagnetization upon each vertex, whereby the production of monopole/anti-monopole pairs locally yield a net magnetization of zero upon the relevant sub-lattice. A second possibility is the production of stripes of alternating magnetization direction, yielding complete demagnetization upon a given sub-lattice. Here magnetic charges can only be found at the stripe ends. A key quantity which will be important in determining the means of demagnetization is that of the monopole effective chemical potential, which quantifies the extent to which monopoles remain closely correlated. This is defined as \({\mu }^{\ast }=\mu /u\), where \(u={\mu }_{0}{Q}^{2}/4\pi a\). We note that when this value approaches half the Madelung constant (for diamond lattice, M/2 = 0.819)26, a highly correlated monopole crystal is energetically favourable and hence is a possible state during the field-driven dynamics. Within a simple dipolar model, for a single spin-flip upon a coordination four vertex (Fig. 6a) we calculate (see ‘Methods’) \({\mu }^{\ast }=1.03\). Surface energetics restrict magnetic charges upon coordination two vertices, so we must consider both the high-energy, coordination two intermediate state, modulated by the factor \(\alpha\), and the final state in which some energy has been spent separating the monopoles from this intermediate state (Fig. 6b). We note that due to the micromagnetic nature of these nanowires, which switch via domain wall motion, this intermediate state must be surpassed.

Fig. 6: Energy landscape for magnetic charges upon the 3DASI lattice.
figure 6

a Energy landscape for monopole nucleation upon the L2 sub-lattice. A single nucleation event from an ice-rule state (a), yields a higher energy state (b) whereby monopoles of opposite magnetic charge reside upon adjacent coordination four vertices. b Energy landscape for monopole nucleation upon L1. Here, a single nucleation event from an ice-rule state (a), yields an intermediate high-energy state (b), which is not observed in the experimental data. A further spin-flip then yields monopoles of opposite magnetic charge, separated by two wires (c), overall reducing their magnetic Coulomb interaction.

Full size image

The intermediate state requires an increase in energy per monopole of \({E}_{{{\mathrm{int}}}}^{\ast }=5.16\). Though the system must clear this energy barrier to transition to a more favourable state, it is more conventional to only consider the chemical potential with respect to the final state. Due to the less favourable Coulomb interaction of the monopoles, the effective chemical potential to produce a monopole across an L1 coordination two vertex is \({\mu }^{\ast }=1.22\), overall yielding a larger fraction of uncorrelated charges. It is therefore this higher \({\mu }^{\ast }\) and the presence of a significant energy barrier for formation of magnetic charge upon a coordination two vertex that overall yields a large proportion of uncorrelated monopoles upon L1.

A key question that remains is the magnitude of surface energetic factor (α) and why such large values are required in MC simulations (α = 6.4) when compared to the magnitude implied by micromagnetics. The surface energetics in these systems arises due to a difference in how the magnetic charge is distributed for two-way and four-way junctions24. In both cases, this will be dictated by a balance between exchange and dipolar energies. For coordination two vertices, the magnetic charge density has a smaller area over which it can distribute, since it can only travel along two wires, producing an unfavourably large energy for monopoles upon the vertex. In contrast, the coordination four system allows the magnetic charge to spread across the vertex area into all four wires, reducing the energy and yielding a stable monopole configuration24. It is important to note that even when α = 3.2 (value indicated by micromagnetic simulations) the resulting MC simulations still bear a far closer resemblance to experiments than when enhanced surface energetics are not considered (α = 1), in terms of string length, monopole density, and density of charges upon surface coordination two vertices. However, increasing α beyond the value predicted by micromagnetic simulations yields an even closer resemblance to experiments, due to fundamental differences in the two methods. Specifically, the MC simulations use a compass needle model, where the magnetic charge associated with each wire is distributed evenly across each needle, effectively reducing the energy barrier for surface charges to form. Therefore, a greater value of α is required to suppress surface charges and hence approximate the experimental observations.

In conclusion, we have demonstrated the fabrication of a 3DASI system, where the magnetic configuration upon the upper two nanowire layers can be determined. We find a striking difference in the field-driven magnetic monopole transport along two principal axes. With a field applied along the projection of surface termination, magnetic imaging shows a low number of uncorrelated monopoles during the switching, which are always found at coordination four vertices. Applying a field along the projection of L2 yields large numbers of correlated monopoles. Micromagnetic and Monte–Carlo simulations, supported by simple calculations within a dipolar framework, suggest it is the difference in effective chemical potential, as well as the energy landscape experienced during surface monopole dynamics, which accounts for the measured differences. We anticipate that our study will inspire a new generation in artificial spin-ice study whereby the ground state in these 3DASI systems is explored as a function of key parameters such as magnetic moment and lattice spacing. Ultimately, this may also yield the realisation of monopole crystals as predicted in bulk spin-ice26 or bespoke spin-ice ground states only possible in artificial systems of novel 3D geometry. By utilizing a full suite of magnetic imaging techniques including MFM, nanoscale ballistic sensing27 and novel synchrotron-based methods28, it is hoped that full 3D characterization of the bulk and surface will soon be possible.

Methods

Fabrication

Diamond-bond lattice structures were fabricated upon glass coverslips via two-photon lithography (TPL). Substrates are first cleaned in acetone, followed by isopropyl alcohol (IPA), and dried with a compressed air gun. Next, droplets of Immersol 518 F immersion oil and IPL-780 photoresist are applied to the lower and upper substrate surfaces, respectively. Using a Nanoscribe TPL system, a polymer scaffold in the diamond-bond lattice geometry is defined within the negative tone photoresist, of dimensions 50 µm × 50 µm × 10 µm. Samples are developed in propylene glycol monomethyl ether acetate for 20 min, then 2 min in IPA, to remove any unexposed photoresist. Once again, the samples are dried with a compressed air gun.

Using a thermal evaporator, a uniform 50-nm film of Permalloy (Ni81Fe19) was deposited on the samples from above, yielding a magnetic nanowire lattice upon the polymer scaffold. This deposition requires a 0.06 g ribbon of Ni81Fe19, washed in IPA, and evaporated in an alumina-coated molybdenum boat. A base pressure of 10−6 mBar is achieved prior to evaporation, the deposition rate is 0.2 nm/s, as measured by a crystal quartz monitor.

Scanning electron microscopy

Imaging was performed using a Hitachi SU8230 SEM with an accelerating voltage of 10 kV. Images were taken from top view as well as at a 45° tilt with respect to the substrate plane.

Magnetic force microscopy

MFM measurements were performed in tapping mode using a Bruker Dimension 3100 Atomic Force Microscope. Commercial low moment MFM tips were magnetised along the tip axis with a 0.5 T permanent magnet. Once mounted, uniform magnetic fields could be applied parallel to each sub-lattice using a bespoke quadrupole electromagnet, which was fixed upon the surface of the AFM stage. During the application of a field, the MFM tip was positioned several mm above the scanning height, such that the tip magnetisation was not influenced. MFM data was taken at a lift height of 100 nm. Prior to capturing MFM images, feedback settings were carefully optimised to ensure sample topography was being accurately measured on the three uppermost lattice layers (L1, L2, L3).

In order to probe the transport of magnetic charge upon the 3DASI surface, the system was placed into a well-defined state by saturating the array along a principal direction (Hsat obtained via optical magnetometry, see Supplementary Fig. 3). MFM images were then obtained after successive field increments in the reverse direction. MFM measures the stray field gradient d2Hz/dz2 due to magnetic charges and hence is an ideal methodology to visualise such transport across the surface25.

Finite element simulations

Micromagnetic simulations of bipod and tetrapod structures were carried out with NMAG28, using finite element method discretisation. Geometries possessing wires with a crescent-shaped cross-section were designed such that the arcs subtend a 160° angle. The inner arc is defined from a circle with 80-nm radius corresponding to the 160 nm lateral feature size of the TPL system. The outer arc is based on an ellipse with an 80 nm minor radius and 130 nm major radius, yielding a thickness gradient with a peak of 50 nm. The length of all wires is set to 780 nm, due to computational restraints. All geometries were meshed using adaptive mesh spacing with a lower limit of 3 nm and upper limit of 5 nm. Simulations numerically integrated the Landau–Liftshitz equation upon a finite element mesh. Typical Ni81Fe19 parameters were used, i.e. MS = 0.86 × 106 Am−1, A = 13 × 10−12 Jm−1 with zero magnetocrystalline anisotropy.

Monte–Carlo simulations

Each nanomagnet is modelled as an infinitesimally thin compass needle with a uniform, linear magnetic moment density mL, with exceptions to this being made at coordination number two vertices. The moment orients along the long axis of the island. The interaction between compass needles is equivalent to two equal and opposite magnetic charges with charges m/L placed at their ends with exceptions for the coordination of two vertex energy. The coordination four energy calculated by the micromagnetic simulations corresponds to compass needles with L = 0.92a while the micromagnetic energies for coordination two imply an interaction strength 3.2 times stronger than those between other charges. The compass needle energy obeys Coulomb’s law with corrections for experimental considerations:

$${E}_{{ij}}={\alpha }_{{ij}}{d}_{i}\frac{{\mu }_{0}{{\rm{m}}}^{2}}{4\pi {L}^{2}}\left[\frac{1}{\left|{{\boldsymbol{r}}}_{{ai}}-{{\boldsymbol{r}}}_{{aj}}\right|}-\frac{1}{\left|{{\boldsymbol{r}}}_{{ai}}-{{\boldsymbol{r}}}_{{bj}}\right|}-\frac{1}{\left|{{\boldsymbol{r}}}_{{bi}}-{{\boldsymbol{r}}}_{{aj}}\right|}+\frac{1}{\left|{{\boldsymbol{r}}}_{{bi}}-{{\boldsymbol{r}}}_{{bj}}\right|}\right]-{{d}_{i}{\boldsymbol{m}}}_{i}\cdot {\boldsymbol{B}}$$
(1)

Where \({{\boldsymbol{r}}}_{{ai}}\) and \({{\boldsymbol{r}}}_{{bi}}\) are the locations of the positive and negative magnetic charge on the ith nanomagnet, \({\boldsymbol{B}}\) is the external field applied either in the L1 or L2 direction with simulated magnitudes ranging between 1.89 and 3.77 mT, \({{\boldsymbol{m}}}_{i}\) is each nanomagnet’s magnetic moment with amplitude m = MV (with \(M\) being the saturation magnetization, and V the nanomagnet volume), \({\mu }_{0}\) is the magnetic permeability, and L = 1000 nm is the island length. The magnetization was chosen to be \({\boldsymbol{M}}=850\) kA m−1 in agreement with previous studies on nanoislands fabricated in the same manner. \({\alpha }_{{ij}}\) is a factor which increases the interaction strength between \({ij}\) pairs at coordination number 2 vertices on the L1 sub-lattice with respect to \({ij}\) pairs at all other vertices, this captures the enhanced surface energetics indicated by micromagnetic simulations. Reversals were simulated with \({\alpha }_{{ij}}\) = 1, 3.23 and 6.45 at coordination number 2 vertices on L1, \({\alpha }_{{ij}}\) = 6.45 was found to yield the closest agreement to experimental observations, all other vertices were consistently defined with \({\alpha }_{{ij}}\) = 1. Site disorder \({d}_{i}\) is drawn from the distribution \(P(d)=\frac{1}{\sigma \sqrt 2\pi } e^{-\frac{(d-1)^{2}}{2\sigma^{2}}}\) (\(\sigma =30\)% is found to yield good agreement with experimental data in this case). This disorder arises due to subtle variations in 3D nanowire geometry. Systems of the same dimensions as the experiment, one unit cell deep and 15 by 15 unit cells wide and long, were simulated with 20 replicas apiece.

The simulations began in the saturated state seen in our experiments. A random spin was selected and the energy to flip that spin was calculated. The corresponding spin-flip was carried out if the energy lost exceeded a threshold energy corresponding to the magnetic coercivity, an algorithm equivalent to the Metropolis method with zero temperature and used in prior spin-ice studies17,29. Flips were attempted 10 times the number of total spins, sufficient for equilibration. The resulting arrow maps were plotted in Fig. 5a–f and Supplementary Fig. 12. In addition, the number of excited states was recorded from these final arrow maps and plotted as a function of field in Fig. 5g.

Dipolar approximation calculations

One can define an effective chemical potential \({\mu }^{\ast }=\mu /u\), where \(u={\mu }_{0}{Q}^{2}/4\pi a\) and \(Q=2\frac{m}{a}\), the smallest non-zero quanta of magnetic charge attainable in our system. Here \(\mu\) is calculated within the dipolar approximation. Magnetic moments are located upon a diamond-bond lattice. The energy of interaction between moments can then be approximated as:

$${E}_{12}=\frac{u}{4}\frac{\left|{\widehat{{\boldsymbol{m}}}}_{1}\cdot {\widehat{{\boldsymbol{m}}}}_{2}-3\left({\widehat{{\boldsymbol{m}}}}_{1}\cdot \widehat{{\boldsymbol{r}}}\right)\left({\widehat{{\boldsymbol{m}}}}_{2}\cdot \widehat{{\boldsymbol{r}}}\right)\right|}{{\left|\frac{r}{a}\right|}^{3}}$$
(2)

This choice of deunitisation places the physical value of the lattice parameter a in u. All geometric factors are then derived from a lattice where the lattice constant is one. The energy landscape for creation of an L2 monopole is shown in Fig. 6a. Considering any two spins within a coordination four geometry yields \({E}_{12}=\frac{5}{4}\sqrt{\frac{3}{8}}u.\) For an ice-rule state there are four low energy pairs and two high-energy pairs yielding \({E}_{{{\mathrm{ice}}}}=-\frac{5}{2}\sqrt{\frac{3}{8}}u\) (point a in Fig. 6a). A doubly charged monopole can be created by flipping a single spin. Each monopole state has three low energy pairs and three high-energy pairs, making \({E}_{{{\mathrm{monopole}}}}=0\) (point b in Fig. 6a).

It should be noted that this discussion of chemical potential typically ignores the next-nearest neighbour spin interactions, or, equivalently, the Coulomb interactions between the generated charges. Monopoles in a Coulomb phase dissociate in a rare gas, leading to an average interaction energy of zero. However, a Coulomb gas is not the final product in this system, so we must consider the pair of charges to begin a lattice spacing apart and interact as dictated by the dumbbell model. That is, long-range interactions take the form \({E}_{{lr}}=\frac{u}{{r}_{{\mathrm{charge}}}/a}\), where \({r}_{{{\mathrm{charge}}}}\) is the distance between the monopoles. At a single lattice spacing, this reduces the energy of charge pair formation by \(u\), or \(\frac{u}{2}\) for each charge. The chemical potential can thus be obtained via \(\mu =({E}_{{\mathrm{monopole}}}-\frac{u}{2}-{E}_{{\mathrm{ice}}})=(\frac{5}{2}\sqrt{\frac{3}{8}}-\frac{1}{2})u\). Therefore, the effective chemical potential, \({\mu }^{\ast }\), for a coordination four geometry is \((\frac{5}{2}\sqrt{\frac{3}{8}}-\frac{1}{2})\) or about 1.0309. In units of Kelvin this yields ~600 K, as expected for a system in the frozen regime.

Now turning to the creation of a monopole upon the L1 sub-lattice, the overall energy landscape is shown in Fig. 6b. Surface energetics restrict the presence of magnetic charges on the coordination of two vertices. Imposing this constraint now requires the resulting monopole–antimonopole pair to separate now \(\frac{2\sqrt{2}}{\sqrt{3}}\) lattice spacings apart. Using the dumbbell model for long-range interactions, the energy required to produce a surface monopole rises by the difference between the interaction of monopoles a lattice spacing apart and monopoles separated further. This instead only drops the cost of monopole formation by \(u\frac{\sqrt{3}}{2\sqrt{2}}\) for two charges, yielding an effective chemical potential for each charge, \({\mu }^{\ast }\), of \(\frac{5}{2}\sqrt{\frac{3}{8}}-\frac{\sqrt{3}}{4\sqrt{2}}\), or ~1.2247 (direct transition between point a to point c in Fig. 6b). The chemical potential calculation upon the L1 sub-lattice does not take into account the detailed dynamics of how these micromagnetic wires switch. Since switching is dictated by domain wall propagation, an intermediate high energy coordination two monopole, though not necessarily stable, needs to be surpassed as part of the monopole nucleation dynamics. The emergence of a coordination two monopole is from a single favourable interaction, −\(\frac{5}{4}\sqrt{\frac{3}{8}}\alpha u\), to a single unfavourable interaction, \(\frac{5}{4}\sqrt{\frac{3}{8}}\alpha u\), each of which is modified by the α factor described in the Monte–Carlo simulations. Again, including the Coulomb correction, this requires an energy change per monopole of \(\left(\frac{5}{2}\sqrt{\frac{3}{8}}\alpha -\frac{1}{2}\right)u\) on the bipod and \(\left(\frac{5}{2}\sqrt{\frac{3}{8}}-\frac{1}{2}\right)u\) on the tetrapod. The average of the two is the intermediate energy per monopole, \({E}_{{{\mathrm{int}}}}=\left(5\frac{\alpha +1}{4}\sqrt{\frac{3}{8}}-\frac{1}{2}\right)u\) (point b on Fig. 6b). For our experimentally confirmed \(\alpha =6.4\), this implies a transition state that costs \({E}_{{{\mathrm{int}}}}^{\ast }=5.1644\).

Magneto-optical Kerr effect magnetometry

A 0.5 mW, 637-nm wavelength laser was expanded to a diameter of 1 cm, passed through a Glan–Taylor polarizer to obtain an s-polarized beam, then focused onto the sample using an achromatic doublet (f = 10 cm), to obtain a spot size of ~10 μm2. During the source-to-sample path the laser is attenuated, approximately reducing the power by a factor of 4. The reflected beam was also collected using an achromatic doublet (f = 10 cm) and passed through a second Glan–Taylor polarizer, from which the transmitted signal was directed onto an amplified Si photodetector, yielding the Kerr signal. After magneto-optical Kerr effect data was captured from the nanowire lattice, a second dataset was obtained from the substrate film. The film data was scaled to the lattice data and subtracted off. This corrected for any Kerr signal contributions originating from the film, during the lattice measurements.