Exploring the phases of 3D artificial spin ice: From Coulomb phase to magnetic monopole crystal

Artificial spin-ices consist of lithographic arrays of single-domain magnetic nanowires organised into frustrated lattices. These geometries are usually two-dimensional, allowing a direct exploration of physics associated with frustration, topology and emergence. Recently, three-dimensional geometries have been realised, in which transport of emergent monopoles can be directly visualised upon the surface. Here we carry out an exploration of the three-dimensional artificial spin-ice phase diagram, whereby dipoles are placed within a diamond-bond lattice geometry. We find a rich phase diagram, consisting of a double-charged monopole crystal, a single-charged monopole crystal and conventional spin-ice with pinch points associated with a Coulomb phase. In our experimental demagnetised systems, broken symmetry forces formation of ferromagnetic stripes upon the surface, a configuration that forbids the formation of the lower energy double-charged monopole crystal. Instead, we observe crystallites of single magnetic charge, superimposed upon an ice background. The crystallites are found to form due to the intricate distribution of magnetic charge around a three-dimensional nanostructured vertex, which locally favours monopole formation. Our work suggests that engineered surface energetics can be used to tune the ground state of experimental three-dimensional ASI systems.

The many-body interaction of dipoles is crucial to understanding a diverse range of phenomena across physics, with its long-range anisotropic nature yielding a wealth of fascinating phenomena.For example, dipolar interactions can yield novel vortex stripes in an ultracold quantum gas 1 , a low temperature residual entropy in frustrated condensed matter systems 2 and Rosensweig instabilities in ferrofluids 3 , yielding self-organised surface structures.The pioneering work of Luttinger and Tisza 4 provided a foundation for understanding dipolar ordering in simple lattice geometries, but this was only extended recently to arbitrary geometries 5 .To date, the experimental placement of dipoles into complex 3D arrangements has been lacking, with scientists mainly relying upon arrangements provided by condensed matter systems.One such model system, known as spin-ice 6 , has been studied intensively and has allowed systematic study of frustration and associated emergence 7 .These systems consist of rare earth moments on corner sharing tetrahedra.The Hamiltonian consists of dipolar and exchange terms and the ~10 Bohr magneton moment means that dipolar interactions are important in determining the nature of the ground state.Since all pairwise interactions within a single tetrahedra cannot be simultaneously satisfied, the system is geometrically frustrated.This yields a local ordering principle known as the ice-rule, in which two spins point into the centre of a tetrahedron and two spins point out, yielding a macroscopically degenerate ground state and a residual entropy measured at low temperature.Interestingly, Monte-Carlo (MC) simulations which encompass sufficient dynamics via a loop algorithm find an ordered phase in spin-ice at very low temperatures, which consists of stripes of anti-parallel spins 8 but so far this has not been measured experimentally.
A new framework to understand the physics of spin-ice was later proposed which treated each spin as a dimer consisting of equal and opposite charges 9,10 .Within this framework local excitations above the ice-manifold are magnetic monopoles in the vector fields M and H, since once the chemical potential has been surpassed, they interact via a magnetic equivalent of Coulomb's law.Subsequent studies provided experimental evidence of magnetic charge transport in bulk spin-ice materials 11,12 .The ground state of spin-ice and the associated dynamic route can then be considered within the framework of magnetic charge, where the ratio of the chemical potential to the magnetic Coulomb energy of a nucleated pair is an important quantity 13 .When this effective chemical potential approaches a value of half the Madelung constant (M/2 = 0.819 for a diamond lattice), a magnetic charge crystal is expected, whereby charges of alternating polarity are tiled throughout the structure 13 .For the canonical spin-ice materials, the effective chemical potential is 1.42, suggesting the monopoles are free to propagate through the system, yielding a disordered spin-ice phase.To observe charge-ordered states in bulk solid-state systems, one needs to find systems with specific material properties.One example, the spin-ice candidate Nd2Zr2O7 has recently shown charge crystal behaviour, combined with disordered spin background, a signature of magnetic fragmentation whereby the local magnetic moment splits into divergence-full and divergence-free parts 14 .A tuneable, engineered system has the capability to explore this phase space systematically.
Artificial spin-ice materials are arrays of lithographically patterned single-domain nanomagnets 15,16 .As such they are a powerful means to explore ordering in dipolar systems by design.Initial studies focussed upon simple square 15 and Kagome arrays 17 , which has subsequently been extended to a wide range of 2D geometries providing a means to explore a variety of model spin systems in statistical physics and more exotic phenomenon such as topological frustration in the Shatki lattice 18 and superferromagnetism in pinwheel lattices 19 .To date, most ASI studies have focussed upon 2D systems due to ease of fabrication but interest has spanned into layered systems 20,21 with both theoretical and experimental studies investigating how these can be used to realise a range of ground states including model vertex systems 22 and superlattice structures 23 .
The advent of three-dimensional lithography now allows the creation of lattices that directly mimic bulk spin-ice geometries [24][25][26] , but with tunability to control factors such as magnetic moment and lattice spacing.Such 3D artificial spin-ice (3DASI) systems, which have a Hamiltonian governed purely by dipolar energetics have only recently been experimentally realised and through simple linear field driving protocols, magnetic charge has been directly observed across the 3DASI surface 24 .
In this article we first use finite temperature MC simulations to carry out a detailed mapping of ordering in idealised 3DASI systems within a diamond-bond lattice geometry.We find a rich phase diagram consisting of a double-charged monopole crystal, single-charged monopole crystal and a spinice phase.We move on to measure the demagnetised state in an experimental 3DASI system and find evidence of an out-of-equilibrium state, whereby crystallites of magnetic charge are superimposed upon an ice background.

Simulating the phase diagram of an idealised 3D artificial spin-ice
Figure 1a and 1b shows a schematic of the simulated unit-cell geometry.Compass needle dipoles are placed upon a diamond-bond lattice, which has a lateral extent of 15 x 15 unit cells and a thickness of a single unit cell.To aid in discussion, we define a series of sub-lattices which are labelled L1 -L4.The upper surface terminates in coordination-two vertices (L1), below which two layers of coordinationfour vertices are found (L2, L3).Finally, the lower lattice surface again terminates in coordination-two vertices (L4).This geometry matches our experimental 3DASI system.The compass needle model (see Methods), is equivalent to treating magnetic dipoles as two magnetic monopoles with a small variable separation.We use a metropolis algorithm to determine the ground state of the system as a function of the dipole length (b), with a fixed lattice spacing (a=1), over a range of temperatures.and q=[-1,-1].At b=1, low temperature, the ground state entropy s0 of spin ice is evident (Fig 1f).In the Methods we calculate s0 analytically using two models: first, using Pauling's method of independent tetrahedra which is well tested in bulk spin ice.Second, by assuming that the surfaces order first and constrain subsequent layers.Fig 1f shows a closer agreement with the latter model, a fundamental difference between the bulk and slab geometries.As b decreases, the frustration and ground state entropy disappear.
Reducing b lowers the chemical potential and in the low temperature regime this yields a phase transition to a double charge crystal (CII).Of particular interest is how such a crystal forms whilst constrained to an odd number of charge layers.The state is characterised by ±2q charges upon surface coordination-two vertices (L1 and L4) and ±4q charges upon L1/L2 and L3/L4 coordination-four vertices, as portrayed in Exploring the ordering in experimental 3D artificial spin-ice systems A 3DASI system was fabricated to explore the extent to which the idealised theoretical phase diagram can be captured experimentally.The system was fabricated using a combination of two-photon lithography and evaporation (See Methods) 24,28 .Figure 2a shows a scanning electron microscopy (SEM) image of the array which takes a diamond bond lattice geometry and has a lateral extent of 50 µm x 50 µm.Figure 2b shows a zoomed top-view, false-colour SEM image with the upper four sublattices labelled (L1 -L4).As in the simulated systems, the lattice terminates in coordination-two vertices upon the surface, with typical coordination four vertices found below at the L1/L2 and L2/L3 junctions.The lower L4 sub-lattice, again terminates in coordination-two vertices.
Our previous work has shown that individual nanowires are single domain and magnetic force microscopy (MFM) can be used to determine the contrast for different vertex types 24 .We now exploit this to determine the demagnetised configuration obtained in 3DASI systems.Note, due to the limited resolution of MFM with lift height, we are only able to measure contrast upon the upper three layers, L1-L3.MFM was performed over large portions of the lattice after planar demagnetisation protocols (See Methods).All vertex types observed in previous experiments 24 , including ice-rule vertices with zero magnetic charge and monopole states with magnetic charge Q=±2q are again observed (Fig 2d).The demagnetised array also contains previously unseen monopole states of charge Q=±4q, as can be seen in  , shows that this is the case for the entire measured area, with coordination-two monopoles being very rare and only observed upon <1% of vertices consistent with previous work 24 .Over large regions of the measured area (~20 %), including in the charge crystallite regions, the L2 sub-lattice is found to host anti-ferromagnetic ordering (Extended Data Fig 2a,b).Breaks in the antiferromagnetic ordering upon L2, via short ferromagnetic strings occurs frequently, with frequency decaying with string length (Extended Data Fig 2c).Interestingly, we find that breaks in the antiferromagnetic ordering often occurs to mitigate the formation of +/-4q charges.We note that since the configuration of the charge crystallites observed experimentally (CIE) has ferromagnetic stripes on L1, it is distinct to the CI charge crystal seen in simulations.
Between areas of magnetic charge crystallite, large patches of the ice phase are observed, as shown by the orange region in

Magnetic charge crystallite formation
We now discuss the observed experimental configuration in terms of the states predicted by MC simulations.For the real experimental systems studied here, the scaled needle length (b) depends upon the vertex type (Extended Data Fig 7), due to the presence of domain walls close to the vertex.When considering all ice-rule vertices, an average b of 0.89 is obtained, suggesting a Q=±4q monopole crystal would be expected as the ground state.However, a set of Q=±2q crystallites form, superimposed upon an ice background.A number of factors may account for this discrepancy.Previous work has suggested that in experimental 3DASI systems, magnetic charges upon surface coordination two vertices are very unfavourable with micromagnetic calculations of single vertices suggesting such excitations cost a factor of three larger than coordination four monopoles 24 .The immediate implication of this is that ferromagnetic stripes upon L1 will forbid the formation of a Q=±4q monopole crystal, apart from regions with local disorder.This is reinforced by the deterministic demagnetisation protocol which favours the formation of ferromagnetic stripes upon the surface.Given this constraint, the system can only form a single charge crystal.However, the formation of charge crystallites via a demagnetisation routine remains surprising and has not been seen previously in either pristine, traditional 2DASI or more exotic layered 2.5D systems.In the former case, charge crystals can be formed in modified square ASI by utilising an MFM tip to selectively switch islands 29 but demagnetisation of conventional square systems yields a low magnetisation, disordered ice phase with low frequencies of monopole excitations 15 .For pristine Kagome systems, demagnetisation yields a 2-in/1-out ice-rule throughout the lattice 17 with only thermally annealed systems yielding some degree of charge ordering 30 .Modifications of the Kagome geometry, either by tuning island lengths within a single unit cell 31 , or by placing exotic nano-bridges at vertices, can also yield charge ordering 32 .
Considering the dynamics of the demagnetising protocol and starting in saturating fields, the system becomes uniformly tiled in type II vertices.Though these are the lowest energy state for single vertices 24 , the net magnetisation makes these less favourable globally.The effective chemical potential upon L1 24 as modified by surface energetics (µ*=1.22)means that deconfined monopoles nucleate and propagate for each 180 degree rotation of the field.At threshold fields, nucleation events upon L1 become less likely, leaving long ferromagnetic strings as observed in the experimental data.The effective chemical potential upon the L2 sub-lattice 24 , within a simple dipolar approximation is lower (µ*=1.03)and favours the local production of correlated charge pairs (type III vertices).It is interesting to note that type III vertices also have a slightly lower b value (0.78) due to the complex 3D distribution of magnetic charge around the vertex.When taking into account the reduced separation between nucleated charges, this reduces the effective chemical potential and yields (See methods) a value of µ*=0.91,approaching the critical value of M/2=0.819.The implication of this is that it is very favourable for a monopole pair to nucleate and remain correlated.Once a single pair is formed, local vertex-vertex Coulomb interactions are minimised by tiling charges of opposing sign.The residual ice-rule regions reflect regions which have not yet equilibrated.It is possible that longer or more complex 3D demagnetisation protocols will promote more efficient exploration of the energy landscape, allowing such ice regions to be further minimised.Altogether, returning artificial spin ice to its three-dimensional origins unlocks previously inaccessible exploration of phase space.Not only is the vertex symmetry restored, but the charge excitations are brought close enough together to encourage never before seen charge crystallite formation from a charge neutral background.We anticipate that fine control of 3D vertex geometry and NiFe thickness will allow suppression of surface energetics and together with an exploration of more complex demagnetisation protocols, or thermal relaxation will allow a realisation of the double charged crystal.It is also expected that more sophisticated synchrotron techniques 35 may also allow imaging of systems greater than one-unit cell in thickness.
• Let S be the total entropy of the system, and s=S/(5N).
Paramagnet: we use the high-temperature paramagnetic phase to constrain the entropy in our numerical calculations.Each domain's orientation is independent.In layers 1 and 5 there are two edges per vertex, giving s1,5 = log 2 2 , and in the other layers there are 4 edges per vertex giving s2,3,4 = log 2 4 .Therefore,  = (16/5)log (2) ≈ 2. 22  (1) Spin ice: for the ice state without field annealing (the physically relevant case around b/a ≈ 1) we have a net charge of zero at every vertex.Each L1 vertex therefore has a net magnetization (1-in 1-out means both domains align).It appears the L1 vertices are completely uncorrelated, even along a single L1 line.Therefore, there are 2 choices per vertex, and s1 = s5 = (1/5) log (2) .Each L2 vertex now has two of its domains fixed by L1 vertices, giving no freedom in these two domains.There are two remaining independent choices per vertex, giving s2,4 = (1/5)log (2) .L3 is then completely constrained by L2 and L4 .Therefore S3 =0 .Overall,  = (4/5)log (2) ≈ 0.55 This agrees with our numerically calculated result to within standard error.The value differs from the Pauling estimate in bulk spin ice; this is because surface energetics dominate in our single-unit-cell slabs.

Monte Carlo Simulations
The interaction energy between two artificial nanomagnets may accurately account for their finite size through the compass needle model.That is, the energy,  !" , between magnets  and  is approximated by considering two point charges at the end of each nanomagnet that interact with Coulomb attraction or repulsion: Here  # is the permeability of free space,  is the nanomagnet's magnetic moment, and  is the nanomagnet length. $! is the position of a magnetic charge, the first index  referring to it being positive and the second index  denoting to which magnet it belongs. !" is the surface energy factor, which for data presented in this publication was set to 1. Since nanomagnet length wildly influences energy scales, all computational energies were normalized by their strongest interactions, such that  ; !" =  !" / %$& .
From this energy we can see that increasing length of the magnets increases nearest neighbour dominance.It's worth noting that the exact distribution of charge and, therefore, precisely what the "length" of the magnets is, depends largely on the details of the nanomagnet's geometry and domain wall arrangement.In this study, this energy is used in the evaluation of a metropolis method Monte Carlo analysis.

Effective Chemical Potential Calculations
The chemical potential of a coordination-four vertex, upon a diamond bond geometry has previously been calculated within a dipolar framework.The energy between any pair of dipoles can be written as: Where m represents the magnetic moment unit vector, r is the moment separation, a is the lattice constant and u is the Coulomb energy between charges: with  = 2/ .One can then simply write the chemical potential as the energy difference between a monopole and an ice-rule state, offset by the magnetic Coulomb interaction, between created charges: with Assuming a perfect dipolar model whereby the charges are separated by a single lattice spacing yields a chemical potential  of 1.03u.The effective chemical potential is therefore However, in our real experimental system the charge separation in the monopole state is reduced, with  >@$2A= ≈0.8a, yielding a reduced  * = 0.91.This locally promotes the formation of charge crystallites.

Magnetic Charge Crystal Order Parameter
In charge ordered systems, twofold degenerate patterns emerge as the ground states.To measure similarity to these states, we can calculate a charge crystal order parameter defined as: ∆ ! is a template of +1's, -1's, and 0's representing a ground state.This was used to calculate the order parameter for both Monte Carlo simulations, as shown in Fig 1c and for experiments.

Magnetic Structure Factor
In the canonical spin ice materials, spin-flip neutron scattering 12 provides what is probably the clearest evidence of spin ice behaviour.Neutron scattering probes the magnetic structure factor projected along the direction of neutron propagation.In artificial spin ice there is a similar tradition of calculating the structure factors, although neutron scattering is not used as a probe.Instead, the structure factor can be inferred directly by Fourier transforming the MFM image 16,33 .In this work we calculated the magnetic structure factor for spin configurations modelling those in our real lattices, as well as those generated in our Monte Carlo simulations.We calculated the full 3D structure factors before taking the qZ = 0 slice, suitable for modelling what would be seen when Fourier transforming a surface MFM arrow map.

Mean Field Analysis
Considering the system in the dumbbell model approximation, where  != ±2, ±1, 0 is the value of the charge on the th vertex,  !" is the interaction strength between charges, and  is the chemical potential of a charge.We can calculate the Maxwell-Boltzmann distribution in the mean field approximation and motivate how charge ordering differs from spin ice ground states.Taking the change of variables  != Δ ! !, where Δ ! is a general charge ordered ground state, and introducing a perturbative "field" to this variable ℎ which will later be set to zero, The variable is approximated by deviations from its mean value,  != ⟨⟩ +  ! .The energy gained by a charge ordered state is called the Madelung constant, which can be written as and  = − ∑  !" Δ !Δ " " . Substituting then yields from this we can calculate the partition function of a single variable and, because they are independent,  = ( ' ) C .For a pyrochlore lattice,  != ±2 (Ω = 1), ± 1 (Ω = 4), and 0 (Ω = 6) where Ω is the degeneracy.Substituting  = (− + 2), the expectation value of the charge ordering variable is then obtained self consistently from the partition function: This self-consistency equation is relatively standard for mean field theories.At high values of −, effectively equivalent to low temperatures, ⟨⟩ = ±2.As − becomes closer to zero, these values gradually drop until the system is no longer ordered.This ordering transition may be characterized for small ⟨⟩ through the first term of the Taylor series: This is true when ⟨⟩ = 0 and  ≤ −1, meaning below a critical temperature, the system will transition to a nonzero order parameter, corresponding to a charge crystal.That critical temperature is This agrees with the previous experimental results that found a spin ice ground state in systems with a reduced chemical potential greater than D ( and a lack of discrete transition in this regime.The critical temperature also decreases with chemical potential as previously observed.Also, since as temperature approaches zero, the order parameter approaches 2, a double charged crystal is the anticipated ground state.One can justify this by considering the lower entropy of the doubly charged state.Since the experimental system is limited to single charges on the surface, the maximum order parameter we predict for the charge crystal ground state of the pyrochlore thin film with 5 charge sites is  > I = 1.5.

Fabrication of 3DASI lattices
Three-dimensional artificial spin ice lattices were fabricated using two-photon lithography followed by thermal evaporation of Ni81Fe19.The coverslips were cleaned in acetone in an ultrasonic cleaner and then washed by isopropyl alcohol (IPA), after which samples were gently dried using compressed air.The coverslip was prepared for TPL with a droplet of immersion oil on one side and Nanoscribe negative-tone photoresist (IPL-780) on the reverse side.The coverslip was then loaded into a Nanoscribe TPL apparatus, and a fabrication script created a number of diamond-bond lattice geometries, each with varying power and scan speed settings.The dimensions of each created lattice are 50µm x 50µm x 10 µm.The completed sample was developed in propyl glycol monomethyl ether acetate (PGMEA) and then rinsed in IPA.An air gun was then again used to remove excess IPA.The sample was then subject to a 50 nm Ni81Fe19 evaporation, at a base pressure of 1 x 10 -6 mBar.Approximately 0.06g of evaporated permalloy was used to achieve this thickness based on previous depositions.A crystal quartz monitor (QCM) present during evaporation measured the deposited thickness; this was later confirmed with atomic force microscopy measurements.The resultant structure has a diamond bond geometry polymer scaffold with magnetic material upon the upper surface of nanowires forming a crescent shaped cross-sectional geometry.Due to line-of-sight deposition, the magnetic coating creates a 3DASI lattice which is one unit cell in thickness, as described previously 28 .Individual nanowires are single domain and have a crescent-shaped cross-section with effective width of 200nm and length of 866nm.

Experimental demagnetisation of lattices
We used a demagnetising protocol akin to method 1 in a previous publication 34 with a sample rotating at ~1000 revolutions per minute, with axis perpendicular to substrate plane.This effectively yields a rotating magnetic field in the substrate plane.The magnetic field starts at 0 mT and ramps up to 75 mT at 2.5 Ts -1 where it is held for 1 second.After this, the field ramps down at 2.5 T s /' to −75 mT and is held for a second.The field then oscillates, whilst the magnitude decreases stepwise to zero over a period of five days.

Magnetic force microscopy (MFM)
MFM data was captured using a Bruker (Dimension Icon) scanning probe microscope in tapping mode.Ultra-low moment probes were magnetised along the tip axis using a 0.5T permanent magnet.The samples were placed with the L1 sublattice parallel to the probe cantilever with a 45-degree scan angle to the L1 sublattice.MFM data were captured using a 65nm lift-height.Separate scans with reversed tip magnetisation were performed to verify consistency of the contrast, and separate scans with the sample rotated 180 degrees were performed to control for artefacts in the scans.
Figure1aand 1b shows a schematic of the simulated unit-cell geometry.Compass needle dipoles are placed upon a diamond-bond lattice, which has a lateral extent of 15 x 15 unit cells and a thickness of a single unit cell.To aid in discussion, we define a series of sub-lattices which are labelled L1 -L4.The upper surface terminates in coordination-two vertices (L1), below which two layers of coordinationfour vertices are found (L2, L3).Finally, the lower lattice surface again terminates in coordination-two vertices (L4).This geometry matches our experimental 3DASI system.The compass needle model (see Methods), is equivalent to treating magnetic dipoles as two magnetic monopoles with a small variable separation.We use a metropolis algorithm to determine the ground state of the system as a function of the dipole length (b), with a fixed lattice spacing (a=1), over a range of temperatures.Figure 1c and 1d show an overview of the phase diagram as a function of b over a range of temperatures, whilst Figs 1e and 1f show the specific heat Cv and corresponding entropy per site s for four values of b.To facilitate interpretation, we define an order parameter (Mc, see Methods) which quantifies the extent to which a magnetic charge crystal has formed.For lower temperatures, a high b lattice yields strong local Coulomb interactions upon vertices, forcing charge neutrality and a spin ice ground state as can be seen in Fig 1c and 1d.A representative arrow map of the spin-ice state is shown in Fig 1g.Ice vertices dominate the microstate occurring at frequencies reflecting underlying vertex probabilities (ergodic balance).The surface L1 layer forms short ferromagnetic strings as seen in previous theoretical studies 27 .The magnetic structure factor (Fig 1h) shows pinch points associated with a Coulomb phase and signatures of short-range magnetic strings with diagonal lines seen along q=[1,1]and q=[-1,-1].At b=1, low temperature, the ground state entropy s0 of spin ice is evident (Fig1f).In the Methods we calculate s0 analytically using two models: first, using Pauling's method of independent tetrahedra which is well tested in bulk spin ice.Second, by assuming that the surfaces order first and constrain subsequent layers.Fig1fshows a closer agreement with the latter model, a fundamental difference between the bulk and slab geometries.As b decreases, the frustration and ground state entropy disappear.
Fig 1i and Fig 1d.The order parameter (Mc) of this CII state is found to be greater than 0.8, as seen by the yellow region in Fig 1c.A neutral layer is found in the centre, consisting of type I vertices.Notably, the sheet geometry produces a coarse-grained field that is approximately constant with respect to distance.This makes the inclusion of a neutral spacing layer more negligible.The magnetic structure factor (Fig 1j) shows clear Bragg peaks due to antiferromagnetic order and associated charge ordering.With intermediate values of b, and at higher temperature, one of the coordination-four double-charged sheets "spreads" into the neutral layer, creating two consecutive single charge sheets, as depicted in Fig 1k and a cross-sectional view shown in Fig 1d, right-panel.This state is named CI.This increases the entropy of the system while maintaining a relatively favourable environment for charges.As temperature is further increased, a peak in specific heat (Fig 1e), corresponding increase in entropy (Fig 1f) and decrease in Mc (Fig 1c) indicates a phase transition to a paramagnetic state.Overall, the phase diagram described by MC simulations is also captured analytically with a simple mean field analysis (See Methods).

Figure
Figure3ashows an experimental magnetic charge map of a 30µm x 30µm region of the lattice, determined by MFM.Three distinct phases are measured and can be readily identified in the charge map with detailed configuration shown in Figs3b-d.Magnetic charge crystallites can be seen with ±2q tiling, as highlighted by the green box in Fig 3a.An arrow map of a typical charge crystallite region is shown in Fig 3b, which shows that it arises due to two types of distinct ordering.The L1 sub-lattice that consists of alternating coordination two and coordination four vertices is found to order into ferromagnetic stripes.Analysis of the L1 sub-lattice (Extended Data Fig 1), shows that this is the case for the entire measured area, with coordination-two monopoles being very rare and only observed upon <1% of vertices consistent with previous work24 .Over large regions of the measured area (~20 %), including in the charge crystallite regions, the L2 sub-lattice is found to host anti-ferromagnetic ordering (Extended Data Fig2a,b).Breaks in the antiferromagnetic ordering upon L2, via short ferromagnetic strings occurs frequently, with frequency decaying with string length (Extended Data Fig2c).Interestingly, we find that breaks in the antiferromagnetic ordering often occurs to mitigate the formation of +/-4q charges.We note that since the configuration of the charge crystallites observed experimentally (CIE) has ferromagnetic stripes on L1, it is distinct to the CI charge crystal seen in simulations.
Fig 3a, with full representative arrow map shown in Fig 3c.These ice regions are largely composed of type II vertices, which due to a subtle broken symmetry in 3D geometry, are the lowest energy vertex type according to micro-magnetic simulations 24 .Finally, only very small regions of the double-charge (CII) crystallite are observed as shown by purple region in Fig 3a and associated arrow map in Fig 3d.The full measured region is shown in Extended Data Fig 3a, with associated vertex types shown in Extended Data Fig 3b and vertex charge shown in Extended Data Fig 3c.The vertex statistics show a strong preference for type III vertices (61.2%), followed by type II vertices (29.8%).Both low energy type I vertices and high energy type IV vertices are only observed occasionally at 5.3% and 3.6% respectively.As would be expected, our measurements indicate charge neutrality, within error as shown in Extended Data Fig 3c.Overall, the charge order parameter as calculated for simulations takes a value of 0.31, for this experimental system (See Methods).The magnetic structure factor of the entire measured data is shown in Fig 4a, with sub-sets corresponding to individual sub-lattices shown in Figs 4(b-d).Focussing upon the data for all layers (Fig 4a), the presence of intense Bragg peaks can be seen, superimposed upon weaker diagonal lines along q=[1,1].In order to further interpret this data, we deconvolve the layers.The L1 structure factor (Fig 4b) consists of a peak upon q=[0,0], indicative of ferromagnetic order on the surface.Weaker split peaks about q=[1/2,1/2] come about due to presence of longer period domains upon L1, as demonstrated in Extended Data Fig 4. The L2 structure factor (Fig 4c) shows peaks due to both type II tiling as well as the magnetic charge crystallite regions as demonstrated in Extended Data Fig 5. Finally, the L3 structure factor (Fig 4d) shows a diffuse signal, with weak Bragg peaks superimposed.This is consistent with the full arrow map (Extended Data Fig 3a), which shows a mixture of charge-ordered and ice states upon L3.Further breakdown of the structure factor via layer and region can be found in Extended Data Fig 6.

Figure 1 :
Figure 1: Simulating the phase diagram in a 3D artificial spin-ice.(a) A unit cell of the simulated geometry.Spins are placed onto the bonds of a diamond-bond lattice.Magnetic charges are represented as spheres of different colour according to legend.(b) View of unit-cell along [001] direction.Arrows are coloured according to layer with cream denoting the surface termination (L1) with coordination-two vertices, yellow (L2), brown (L3) denoting coordination-four vertices and dark red denoting lower surface termination (L4) with coordination-two vertices.(c) Phase diagram of 3D artificial spin-ice showing charge crystal order parameter (Mc) as a function of

Figure 2 :
Figure 2: Experimental vertex states and measured magnetic force microscopy contrast.(a) A scanning electron microscope image of the 3D artificial spin-ice lattice, scale bar is 20 µm.(b) Top-view, false-colour SEM image of the 3DASI lattice with the individual sub-lattices labelled.Scale bar is 1 µm.(c) Magnetic force microscopy contrast of the vertex types found upon the L1 coordination two vertices.(d) Contrast for vertex types measured at the L1/L2 vertex.Here, both ice-rule (Type I, Type II) and single charged vertices (Type III) are observed.(e) Contrast for vertex types measured at the L2/L3 vertex, which shows a mixture of ice vertices as well as those with single and double magnetic charge (Type IV).

Figure 3 :
Figure 3: Measuring the experimental demagnetised state of a 3D artificial spin-ice.(a) Global magnetic charge map of the measured sample region.Charged regions are represented by colour as according to legend.The map shows examples of single charge crystallite (green outline), the ice phase (orange outline) and the double charge crystal (purple).(b) More detailed arrow map of the experimental single charge crystallite.It can be seen to consist of ferromagnetic stripes on the surface L1 layer with antiferromagnetic ordering upon L2.(c) Ice phase with type II tiling and (d) double charge crystallite, which only occurs at breaks in the L1 ferromagnetic ordering.Arrows represent the magnetisation of the L1, L2, and L3 sublattices.Here, the distinction between surface

Figure 4 : 3 :Extended Data Figure 4 :Extended Data Figure 5 :Extended Data Figure 6 :
Figure 4: Experimental magnetic structure factors.(a) The magnetic structure factor of all sub-lattices superimposed.Clear Bragg peaks can be seen with periodicity in two dimensions.(b) The magnetic structure factor of the L1 sub-lattice.Peaks can be seen at q=[0,0], corresponding to the ferromagnetic ordering upon the surface.The split peaks about q=[1/2, 1/2], are due to domains on L1 with larger periodicity as demonstrated in Extended Data Fig 4. (c) The structure factor for the L2 sub-lattice.Bragg peaks can be seen, resulting from both type II vertices and charge crystallites (d) The structure factor for the L3 sub-lattice.Bragg peaks are again seen and match the periodicity seen for L2 sub-lattice, superimposed upon a diffuse background.
Table showing scaled dipolar needle length for range of vertices.