Persistent dynamic magnetic state in artificial honeycomb spin ice

Topological magnetic charges, arising due to the non-vanishing magnetic flux on spin ice vertices, serve as the origin of magnetic monopoles that traverse the underlying lattice effortlessly. Unlike spin ice materials of atomic origin, the dynamic state in artificial honeycomb spin ice is conventionally described in terms of finite size domain wall kinetics that require magnetic field or current application. Contrary to this common understanding, here we show that a thermally tunable artificial permalloy honeycomb lattice exhibits a perpetual dynamic state due to self-propelled magnetic charge defect relaxation in the absence of any external tuning agent. Quantitative investigation of magnetic charge defect dynamics using neutron spin echo spectroscopy reveals sub-ns relaxation times that are comparable to the relaxation of monopoles in bulk spin ices. Most importantly, the kinetic process remains unabated at low temperature where thermal fluctuation is negligible. This suggests that dynamic phenomena in honeycomb spin ice are mediated by quasi-particle type entities, also confirmed by dynamic Monte-Carlo simulations that replicate the kinetic behavior. Our research unveils a macroscopic magnetic particle that shares many known traits of quantum particles, namely magnetic monopole and magnon.

cial spin ice can be described by the same quantum mechanical treatment as in spin ice, where the Hamiltonian representing magnetic charge interaction is defined by Pauli matrices. 3,4.A direct comparison of magnetic charge relaxation times in spin ice and artificial magnetic honeycomb lattices can shed light on this fundamental problem.From the Coulombic physics perspective, a magnetic charge at the center of an artificial kagome lattice can be mapped onto the net magnetic flux along the bond direction on the parent honeycomb motif. 20,21The correspondence entails the quantum treatment of magnetic charge interactions, which has been used previously to deduce thermodynamic and magnetic properties of the statistical ensemble 20 .A charge defect's motion between vertices would flip the microscopic moment in the nanoscopic honeycomb element, manifesting a magnonic pattern, thereby altering the moment direction and the net charge on a given vertex, see Fig. 1a-c.The charge defect or magnonic charge, q m , cannot pass through a high integer charge (3Q or −3Q), which serves as a roadblock.At the same time, ±3Q charges have a high energy cost, which makes them unstable and triggers the dynamic process.A thermally tunable honeycomb lattice, made of ultra-small (∼11 nm) nanoscopic permalloy element (see Fig. 1d), 22 can host such dynamic event without the application of any external stimulus.It is known to depict a highly disordered ground state, comprised of both ±Q and ±3Q charges, at low temperature. 23The ultra-small element size imparts thermally tunable energetics due to the modest inter-elemental magnetic dipolar interaction energy of ∼ 45 K.We use neutron spin echo (NSE) measurements to obtain direct evidence of self-propelled relaxation of magnonic charge q m , as it provides accurate quantitative estimation of relaxation properties across a broad dynamic range. 24Furthermore, FIG.1: Self-propelled magnetic charge defect relaxation in an artificial honeycomb lattice made of ultra-small permalloy elements.a, Schematic illustration of magnetic charge relaxation processes between honeycomb vertices.The charge defect (qm) dynamics between nearest neighboring vertices corresponds to the reciprocal wave vector q ∼ 0.058 Å−1 .Magnetic charge relaxation occurs simultaneously in all three elements attached to a vertex.b, The dynamic state due to qm kinetics in a honeycomb element resembles magnon dynamics pattern.c, Time-dependent micromagnetic simulations depict magnetization pattern in a honeycomb element at two different instances, t = 0 and t = 363 ps, of qm propagation.d, Atomic force micrograph of artificial honeycomb lattice, with typical element size of ∼ 11 nm (length)×4 nm (width)×6 nm (thickness).
NSE spectroscopy is carried out on two stacks of permalloy honeycomb samples with varying thicknesses of 6 nm and 8.5 nm that render a comprehensive outlook.
In Fig. 2a we show the color plot of NSE data obtained on a 6 nm thick honeycomb lattice with net neutron spin polarization along +Z direction.No meaningful spectral weight is detected in the spin down neutron polarization (see Fig. S2-S3 in SM), 25 which confirms the magnetic nature of the signal.To ensure that the observed signal is magnetic in origin, NSE spectrometer is modified by removing the flipper before the sample (see Methods for detail).The figure exhibits very bright localized scattering pixels that are identified with q = 0.058 Å−1 and 0.029 Å−1 .These q values correspond to the typical relaxation length associated with the distance between the nearest and the next-nearest vertices, 2π/l and 2π/2l, respectively, with l = 11 nm.This immediately suggests the quantized longitudinal nature of magnetic charge defect dynamics in artificial honeycomb lattices.NSE spectroscopy is a quasi-elastic measurement technique where the relaxation of magnetic specimen is decoded by measuring relative change in scattered neutron's polarization via the change in the phase current at a given Fourier time (related to neutron precession).The localized excitations exhibit pronounced sinusoidal spin echo.Fig. 2b,c show the echo profiles obtained on honeycomb samples at room temperature at 0.02 ns Fourier time where a strong signal-to-background ratio is detected.
Typically, the echo intensity from a single wavelength λ follows a cosine function, given by I(φ) = A cos(φλ) + B, where the phase φ is directly related to the phase current dJ.However, if the neutron has a wavelength span of λ avg ± dλ, then the signals from all different wavelengths add up (see SM).In this case, the total intensity follows the following functional relation, 26 As the Fourier time associated to neutron precession increases, the error bar tends to become larger.Fig. 2d shows the spin echo at 0.26 ns Fourier time (see Fig. S4 in SM for other Fourier times).This trend hints at the fast relaxation of the magnetic charge quasi-particles, which is more prominent at a sub-ns time scale.The sinusoidal nature of spin echo is quite evident at other temperatures as well.In Fig. 2e,f, we show the plot of characteristic spin echo profiles at low temperature in both thicknesses at q = 0.058 Å−1 .Similar echo profiles are observed at q = 0.029 Å−1 .The strong spin echo at low temperature suggests the persistence of significant charge dynamics in the nanoscopic system.The origin of charge relaxation at low temperature of T = 4 K is unlikely to be thermal in nature considering that it costs |E Q -E −Q | (∼ 76 K) to create a charge defect.Unlike other nanostructured magnets where the domain wall motion, typically ascribed to the dynamic property, is detectable in applied magnetic field or current only, the absence of any external tuning parameter at low temperature suggests the non-trivial nature of magnonic charge q m relaxation.In total, measurements were performed at six temperatures of T = 4 K, 20 K, 50 K, 100 K, 200 K and 300 K (see Fig. S5 in SM).The quantitative estimation of relaxation times at different temperatures sheds light on the dynamic state of artificial honeycomb lattice.
The estimation of magnetic charge relaxation time is achieved by analyzing the intermediate scattering function S(q, t)/S(q, 0) at different spin echo Fourier times.We show the plot of normalized intensity as a function of neutron Fourier time at q = 0.058 Å−1 in the 6 nm thick sample in Fig. 3a.The normalization is achieved by dividing the observed oscillation amplitude by the maximum measurable amplitude (see SM), which is a common practice in magnetic systems with unsettling fluctuation to the lowest measurement temperature 26 .The normalized intensity reduces to the background level above the spin echo Fourier time of ∼ 0.5 ns, conforming to the most general signature of relaxation process in NSE measurements 26,27 .The quantitative determination of magnetic charge relaxation time involves exponential fitting of the NSE scattering intensity as S(q, t)/S(q, 0) = Ce −t/τm , where C is a constant 28  anism of charge q m .Four important conclusions can be drawn from the estimated value of the relaxation time: first, the obtained value of τ m = 20 ps at T = 300 K suggests a highly active kinetic state in our permalloy honeycomb lattice, typically not observed in nanostructured materials where the domain wall motion is the primary mechanism behind the dynamic property 29 .Moreover, magnetic field or current application is necessary to trigger domain wall dynamics, which is not the case here: the dynamics state of charge q m is entirely self-propelled.Second, the relaxation rate of q m in honeycomb lattice is somewhat thickness independent.In Fig. 3b, we show S(q, t)/S(q, 0) as a function of neutron Fourier time in the 8.5 nm thick lattice.The estimated τ m = 49 (10) ps is comparable to that found in the 6 nm thick lattice (within 2σ).Importantly, we do not observe a drastic difference between the charge dynamics in the thinner and the thicker lattice, which is a strong indication for its quasi-particle characteristic.
Third, the charge defect q m relaxes at the same high rate at low T despite the absence of thermal fluctuation.As shown in Fig. 3c, the estimated relaxation time as a function of temperature suggests the persistence of charge dynamics at lower temperatures.Given the fact that the charge dynamics prevails in the absence of any external stimuli, the observed T -independence is quite significant.The conclusion is that the system lives in a persistent kinetic state due to the magnetic charge defect motion between parent lattice vertices, which is strong evidence for the quantum nature of magnetic charge quasi-particles.
Lastly, the q m relaxation rate in the 6 nm thick honeycomb lattice is comparable to that found in the bulk spin ice material, ∼ 5 ps. 18,27This suggests a correspondence between the dynamic behavior of magnetic charge defects and the monopolar dynamics in the atomistic spin ice.To our knowledge, this is the first observation of its kind, bridging the gap between the bulk spin ice material and the nanostructured spin ice system.Thus, the quasi-particle characteristic of the monopole can be suitably invoked in the artificial lattice.These conclusions are further supported by theoretical analysis, which we now briefly summarize.
In order to confirm the quantum nature of magnetic charge defect quasi-particles, we modeled the system both fully quantum mechanically and fully classically.Motivated by the simplistic behavior of magnetic charge relaxation, unchanged by temperature or system size, we suspect the system exhibits rapid renormalization group convergence to a minimal effective model of behavior, and as such we assume the system should be modeled effectively by a Hamiltonian with minimal features.To this end, for the quantum modeling we treat the magnetic monopoles at each vertex as a pseudospin 3/2 due to the quartet nature of these sites.These pseudospins are coupled minimally via nearest neighbor flip flop terms, i.e.H nn = ij J h S + i S − j where S + i is the spin 3/2 representation of an SU(2) raising operator for the ith site and the sum is over nearest neighbors.We calculated the average relaxation time as a function of temperature for two cases, a 4-site model solved exactly, and a 200 site model with periodic boundary conditions time evolved via dynamic Monte Carlo. 30, Both cases reproduce the temperature independence of relaxation time for negligible barrier height compared to the large exchange constant, see Fig. 4a (see further detail in the SM).
The main result of the quantum Monte-Carlo simulation is that a magnetic charge defect q m travels through the lattice element effortlessly as a result of negligible energy barrier.A similar observation was made in the case of magnetic monopole dynamics in atomic spin ice. 5 If the barrier height were comparable to the exchange constant J, then the relaxation time would increase significantly at low temperature; however, this is not seen in the experimental data.Similarly, if the kinetic behavior of q m were classical in origin, then the simulation would predict a monotonic temperature dependence of the relaxation time above the barrier height.
For the classical modeling, we have treated each joint of the lattice as a classical spin obeying the Landau-Lifshitz-Gilbert equation with a stochastic field modeling temperature dependence 31 .In this case the magnetic charge relaxation time was found to be temperature independent only up to a threshold temperature T c , above which the relaxation time decreases monotonically with increasing temperature.This threshold occurs when the root-meansquare of the thermal field B th is comparable in magnitude to the effective field term associated with shape anisotropy B an .With an estimated effective anisotropy (mostly shape anisotropy) of |K 1 | ∼ 1 × 10 4 J/m 3 , we estimate the threshold temperature to be on the order of T c ∼ 20 K, see Fig. 4b.This is inconsistent with the experimental observations, leading us to conclude that the charge defect q m dynamics has indeed quantum nature.This is the first time that a comprehensive study of magnetic charge relaxation processes in two-dimensional artificial permalloy honeycomb lattice reveals kinetic events with ∼ps temporal quantization.Previous efforts in determining the charge dynamics in artificial spin ice systems focused on the usage of ferromagnetic resonance spectroscopy and optical measurements such as Brillouin light scattering and optical cavity techniques 32,33 .However, these techniques have intrinsic temporal thresholds of 1 ns or higher, which limits the scope of sub-ns inves- tigation of charge dynamics.Additionally, the samples used for dynamic study were created via electron-beam lithography, which results in large element size.In most cases, samples created in this way are athermal, 9 thus significantly affecting the spatially quantized relaxation process between the vertices of the parent lattice.Our honeycomb sample with small inter-elemental energy, utilized in this study, provided an archetypal platform to extract the charge relaxation properties in the absence of thermal fluctuation.The ultra-small element size imparts a thermally tunable characteristic to the lattice.It is also relevant to mention that the macroscopic size honeycomb specimen, used in this study, tends to develop structural domains that vary in size from 400 nm to few µm.However, the NSE study of charge dynamics only involves the nearest and the next nearest neighboring vertices.Therefore, structural domains are not expected to affect the outcome.The NSE technique, with a dynamic range extending to sub-ns timescales, is the most suitable experimental probe to study such delicate dynamic property.
The synergistic experimental and theoretical investigations have not only revealed the sub-20 ps relaxation time associated to magnonic charge dynamics but also elucidated its quantum mechanical characteristics.The thickness independence of the relaxation time τ gives further support to the quasi-particle nature of the charge q m .Despite the classical characteristic of magnetic charge, directly related to magnetic moment M along the honeycomb element length (L) via q = M/L, the observed quantum mechanical nature of the charge defect dynamics is unprecedented.Unlike magnetic monopole or magnon particles in bulk magnetic materials of atomic origin, q m is expected to possess finite 'macroscopic' size.However, the kinetic behavior of q m manifests similarities with the known properties of magnetic monopoles, which is a strong indication of the particle-type character.
The new finding reported here can be utilized for the exploration of 'magnetricity' in artificial kagome ice, originally envisaged in the spin ice magnet. 6Additionally, the quantum nature of magnonic charge can spur spintronic application via the indirect coupling with electric charge carriers. 34In recent times, experimental efforts are made to develop a new spintronic venue in artificial spin ice system. 35A dynamic system of thermally tunable magnetic honeycomb lattice renders a new research platform in this regard.

Methods Sample Fabrication
Artificial honeycomb lattices were created using hierarchical top-down nanofabrication involving diblock copolymer templates [see details in Supplementary Material (SM) 25 ].An atomic force micrograph of the hon-eycomb lattice sample is shown in Fig. 1d.

Neutron Spin Echo Measurements
We performed NSE measurements on permalloy honeycomb lattice samples using an ultrahigh resolution (2 neV) neutron spectrometer SNS-NSE at beam line BL-15 of the Spallation Neutron Source, Oak Ridge National Laboratory.Time-of-flight experiments were performed using a neutron wavelength range of 3.5 -6.5 Å and neutron Fourier times between 0.06 and 1 ns.NSE is a quasi-elastic technique where the relaxation of a magnetic specimen is decoded by measuring the relative polarization change of the scattered neutron via the change in the phase current at a given Fourier time (related to neutron precession).The experiment was carried out in a modified instrumental configuration of the NSE spectrometer where magnetism in the sample is used as a π flipper to apply 180 o neutron spin inversion, instead of utilizing a flipper before the sample.Additional magnetic coils were installed to enable the XYZ neutron polarization analysis, see the schematic in Fig. S1 in SM.Such modifications ensure that the detected signal is magnetic in origin.NSE measurements were performed on two different parallel stacks of 125 and 117 samples of 20 × 20 mm surface size and 6 nm and 8.5 nm thicknesses, respectively, to obtain good signal-to-background ratio.The sample stacks were loaded in the custom-made sample cell, which was mounted in a close cycle refrigerator with 4 K base temperature.NSE data was collected for 8 hrs on the average at each temperature and neutron Fourier time.
The fabrication of the artificial honeycomb lattice involves the synthesis of a porous hexagonal template on top of a silicon substrate, calibrated reactive ion etching to transfer the hexagonal pattern to the underlying silicon substrate and deposition of permalloy on top of the uniformly rotating substrate in a near-parallel configuration (∼ 3 • ) to achieve the 2D character of the system.The porous hexagonal template fabrication process utilizes diblock copolymer polystyrene (PS)-b-poly-4-vinly pyridine (P4VP) of molecular weight 29K Dalton and volume fractions 70% PS and 30% P4VP, which can self-assemble into hexagonal cylindrical structure of P4VP in the matrix of PS under the right condition.A PS-P4VP copolymer solution of mass fraction 0.6% in toluene was spin-coated on a polished silicon wafer at around 2300 rpm for 30 seconds, followed by solvent vapor annealing at ∼ 25 • C for 12 hours.A mixture of toluene/THF (20/80, volume fraction) was used for the solvent vapor annealing.This process results in the self-assembly of P4VP cylinders in a hexagonal pattern in a PS matrix.Submerging the sample in ethanol for 20 minutes releases the P4VP cylinders from the PS matrix, leaving a hexagonal porous template with an average hole center-to-center distance around 31 nm.The diblock template is used as a mask to transfer the topographical pattern to the underlying silicon substrate using reactive ion etching with CF 4 gas.The top surface of the etched silicon substrate resembles a honeycomb lattice pattern, which is exploited to create a magnetic honeycomb lattice by depositing permalloy (Ni 0.81 Fe 0.19 ) in a near-parallel configuration using E-beam physical vapor deposition.Samples were uniformed rotated to achieve uniformity of the deposition.This allowed evaporated permalloy to coat the top surface of the etched silicon substrate only and producing magnetic honeycomb lattice with a typical element size of about 11 nm (length) × 4 nm (width) and controllable thickness.A typical atomic force micrograph of the honeycomb lattice is shown in Fig. 1d.Samples of thickness 6 nm and 8.5 nm were used in this work.

B. Time-of-flight Neutron Spin Echo instrument
The Neutron Spin Echo (NSE) measurements were conducted on an ultrahigh resolution (2 neV corresponding to ∼ 300 ns Fourier time) neutron spectrometer at beam line BL-15 of the Spallation Neutron Source, Oak Ridge National Laboratory.A 3 Å wavelength bandwidth of neutrons (3.5-6.5 Å) was used in the experiment, and the scattered neutrons were collected with a 30 cm × 30 cm position-sensitive 3 He detector 3.9 m away from the sample.The detector has 32 × 32 X, Y pixels laterally to keep track of the scattering angles of the detected neutrons.It also has 42 time of flight (ToF) channels (tbins) that label the time frame of the detected neutron which encodes the neutron's wavelength.Thus, the data collected by the detector is a 32 × 32 × 42 array.By carefully grouping the X, Y pixels and tbins, one can extract echo signals at different Fourier times t [see Eq. (S2) below for the definition] and q values from a single measurement.In this experiment, echo signals at multiple temperatures and nominal Fourier times (the Fourier time associated with the neutron of the maximum wavelength) were measured.
Unlike in conventional NSE, where a π flipper is installed before or after the sample to flip the neutron's polarization, a magnetic sample itself does a π flip of the neutron's polarization due to the nature of the interactions between magnetic moments and neutron spins.In this experiment, the π flipper was removed and additional magnetic coils were installed for three-dimensional polarization analysis.To obtain good signal-to-background ratio, a stack of 125 (117) samples of about 20×20 mm 2 was loaded in a custom-made aluminum sample container.The sample container was inserted into a close cycle refrigerator with a base temperature of 4 K.A schematic diagram of the NSE instrument is presented in Fig. S1.where S(Q, ω) is the scattering function of the sample and t is the Fourier time, defined as with γ n , m n and λ being the neutron's gyromagnetic ratio, mass and wavelength, and h denotes Planck's constant.J = |B|dl is the magnetic field integral along the neutron's path through the precession coil 1 .The field integrals are designed to be the same in both precession coils before and after the sample such that the precession phase acquired by the neutron in the first coil will be exactly recovered at the end of the second coil, given that the neutron does not change its velocity while interacting with the sample.For quasi-elastic scattering, the neutron gains a net phase at the end of the second coil before it reaches the analyzer.By systematically stepping through an additional phase current (convert to an additional field integral) in either the first or the second coil, a cosine modulation of the sample's intermediate scattering function is realized, which is the typical raw data, echo signal, one would analyze in an NSE experiment.
The echo signal intensity for a given neutron wavelength λ has a cosine function shape given by where φ = dJγ n m n /h and dJ is the phase asymmetry between the two precession coils introduced via the scanning phase current.When signals from neutrons of a certain wavelength span λ avg ± dλ are added up, the total intensity resembles a cosine function modulated by an envelope function With the same color scale, while no meaningful spectral weight is detected in the Z− polarization, in the Z+ polarization case, clear spectral weight is observed in the detector region corresponding to q ∼ 0.058 Å−1 and q ∼ 0.029 Å−1 .The intense signal at the left edge of the detector is the direct beam as the instrument hits its geometrical constraint limit.Fig. 2a is obtained by subtracting the Z− intensities from the Z+ intensities.The color scale on the right shows the detected neutron intensity normalized by the proton charge.With the same color scale, while no meaningful spectral weight is detected in the Z− polarization, in the Z+ polarization case, clear spectral weight is observed in the detector region corresponding to q ∼ 0.058 Å−1 and q ∼ 0.029 Å−1 .The intense signal at the left edge of the detector is the direct beam as the instrument hits its geometrical constraint limit.With the same color scale, while no meaningful spectral weight is detected in the Z− polarization, in the Z+ polarization case, clear spectral weight is observed in the detector region corresponding to q ∼ 0.058 Å−1 and q ∼ 0.029 Å−1 .The intense signal at the left edge of the detector is the direct beam as the instrument hits its geometrical constraint limit.Fig. 2a is obtained by subtracting the Z− intensities from the Z+ intensities.

D. Neutron Spin Echo raw data treatment
The echo intensities were obtained by adding up signals in X, Y detector pixels and then summing over different ToF tbins while correcting the phase shifts among the tbins before the summation process.Grouping of the X, Y detector pixels were guided by the false color detector image (Fig. 2a in the main text) as well as the detector pixel image (Fig. S3) to center around the intense scattering signal.The echo profiles of the 8.5 nm honeycomb lattice shown in Fig. 2c,f of the main text were obtained by summing over detector pixels over the range of X = 8 to 10, Y = 10 to 21 in ToF tbin = 4, with a Fourier time ∼ 0.02 ns.In order to make direct comparisons between the 6 nm honeycomb lattice and the 8.5 nm honeycomb lattice, echo profiles of the 6 nm honeycomb lattice (Fig. 2b,d,e of the main text) with the same Fourier time were obtained by summing over detector pixels over the range of X = 7 to 11, Y = 7 to 24 in ToF tbin = 17.
For the calculation of the intermediate scattering function (Fig. 3a,b of the main text), ToF tbins = 10 to 19, 20 to 29, 30 to 39 were explored.The detector pixel groupings used in the calculation of the intermediate scattering functions at each temperature are identical to that used for the echo profiles.Note that a slightly larger detector area is selected for the data analysis of the 6 nm honeycomb lattice.This is because the detected scattering signal from the 6 nm honeycomb lattice is broader than that from the 8.5 nm honeycomb lattice possibly due to a larger variation of the honeycomb element length scale.

E. Calculation of the intermediate scattering function
The goal of every NSE experiment is to obtain the intermediate scattering function S(q, t)/S(q, 0) as a function of the Fourier time t, which can be calculated from the fitted echo amplitude.For non-magnetic samples, where A denotes the fitted amplitude of the echo signal, U and D denote the spin up (non-π flipping) and spin down (π flipping) measurement of direct scattering without precession.U − D measures the maximum obtainable echo amplitude.For magnetic samples, half of the magnetic scattering intensity M/2 is used for normalization 2 .Thus 6 additional polarization measurements along x, y and z directions were performed at each temperature.M and U, D are calculated as The intermediate function is then determined as which is comparable with 66A/(U − D) in this experiment.In the calculation of the intermediate scattering function at a given Fourier time, the same X, Y pixels and ToF channels selection is used for the determination of both the echo amplitude and the normalization factor.Unlike experiments conducted with other spectrometers, the measured dynamic structure factor S exp (q, ω) is a convolution between the true dynamic structure factor S(q, ω) and the instrument resolution R(q, ω), that is S exp (q, ω) = S(q, ω) * R(q, ω).NSE directly probes the intermediate scattering function, which is the Fourier transform of the dynamics scattering function S(q, t) = cos(ωt)S(q, ω)dω, thus S e xp(q, t) = S(q, t)R(q, t).Therefore, in an NSE experiment, the resolution function can be simply divided out to obtain the true intermediate scattering function.At SNS-NSE the resolution of the instrument and elastic scattering contribution is assessed by measuring a perfect elastic scattering sample, usually mounted in the same container as the sample to investigate, measured over the same q-range, tau-range, and wavelength.For soft matter measurements, solid graphite and Al 2 O 3 as well as TiZr are usually used depending on the scattering angles.For paramagnetic NSE measurements we use Ho 2 Ti 2 O 7 , a well-known classical spin ice material frozen below T = 20 K where it exhibits no dynamics.However, for the measurement performed in this work, it was not possible to use Ho 2 Ti 2 O 7 sample as the resolution since it only scatters and produces reliable echoes at high scattering angles, while the magnetic honeycomb samples scatter in the small-angle regime.Another common practice in quasi-elastic techniques is to measure the same sample at T ≤ 4 K   where presumably all dynamics freezes.However, it is well observable in our measurements that even at T = 4 K fast dynamics of the magnetic charges still exists, which is, in fact, one of the main findings of our research.Therefore, we decided not to reduce the data by elastic resolution, since, as explained, we do not have an elastic scatter with frozen dynamics within the temperature range accessible.In support of our decision is the fact that at SNS-NSE below Fourier times t < 1 ns the elastic resolution is predominantly flat (linear).This means that any reduction by resolution will not affect the relaxation observed in the data but only the intensity scaling on y-axis.To demonstrate the above we have collected two measured Ho 2 Ti 2 O 7 resolutions from previous paramagnetic measurements at T = 4 K and 10 K, for 2 different wavelength 6.5 Å and 8 Å.In Fig. S6, these experimental elastic resolution magnetic data are presented, together with their fits and simulated resolution in the soft-matter regime.One can easily observe the linear behavior in the range of t < 1 ns, which is the predominant range of our measurements.Plots of S(q, t)/S(q, 0) at q = 0.058 Å−1 at intermediate temperatures are presented in Fig. S7.

II. THEORETICAL DETAILS
Here we present a theoretical analysis of the ultra-small artificial magnetic honeycomb lattice, with the intention of elucidating the dominant physics at play behind the experimentally observed temperature independence of the magnetic charge relaxation time.
As discussed in the main text, we suspect the system exhibits rapid renormalization group (RG) convergence to a minimal effective model.This intuition is motivated by the simple, consistent behavior of the magnetic charge relaxation over a broad range of temperatures; a feature that is unchanged with at least a modest change in lattice thickness.By "rapid renormalization group convergence" we specifically mean that small changes to a hypothetical exact theoretical description of the experimental system, such as a change in an interaction term coupling strength or a change in system size, would result in no qualitative change to the physically observable behavior of the model, due to the presence of a well isolated and strongly stable fixed point in RG space.Literature on the structurally equivalent Kagome spin ice model 3 supports the contention that the model exhibits rapid RG convergence to the same high and low temperature phases, with or without the inclusion of long or short range interaction terms.As such we assume the system should be modeled effectively by a Hamiltonian with minimal features, e.g. a Hamiltonian that only captures the magnetic charge excitation energy and the dynamical capacity of these charges to transfer to neighboring sites.
In the following subsections, we consider two minimal effective models, one a classical description and one a quantum mechanical description, and compare the excitation relaxation time for each model to the experimental results.

A. Quantum magnetic charge model
The simplest model for a Kagome spin ice involves spins restricted to specific site dependent directions êi which point along three possible axes separated by 120 • .These axes together form a honeycomb lattice structure akin to the herein studied artificial honeycomb lattices.The spins of the Kagome spin ice are coupled via a nearest-neighbor spin exchange interaction where êi • êi = − 1 2 has been used, and the notation . . . is used to indicate a sum over nearest neighbors.The system is frustrated in the case of ferromagnetic coupling, i.e.J 1 > 0, as can be seen by noting that for this sign of J 1 the  S10) is equivalent to an antiferromagnetic honeycomb Ising model.The Hamiltonian (S10) can be simplified further by introducing the magnetic charge at honeycomb vertex α, Q α = ± i∈α σ i .In terms of these operators, the Hamiltonian is equivalent to, up to a constant, Here the possible values of the magnetic charges are Q α = ±1, ±3, with the ±3 states corresponding to the localized excitations of the model.The model described by H 1 is essentially classical, for the same reason that the original Ising model (without a transverse field) is classical, namely because the model lacks any non-commuting observables.This originates in the fact that, in deriving this Hamiltonian, we have neglected the important off axis components of the spins (S9).Including these off axis terms in the exchange coupling of Eq. (S10) results in nearest neighbor spin flip-flop terms (corresponding to a nearest neighbor hop of Q α s), thus restoring the quantum nature of the model.However, instead of simply reintroducing the spin components, we consider a slightly different starting point in hopes of constructing a minimal quantum mechanical model in terms of the Q α 's.To this end, we note that the quartet nature of Q α suggests that the simplest possible treatment of Q α as a quantum operator would be to treat it as a z component of a pseudospin 3/2 representation of SU(2).With this in mind the natural choice of non-commuting observables would be the ladder operators of the same pseudospin 3/2 representation.With a basic nearest neighbor hopping term, we arrive at the minimal quantum mechanical Hamiltonian where we define q α = 1 2 Q α in order for q α to exhibit the traditional spin 3/2 representation eigenvalues of ± 1 2 and ± 3 2 , and the ladder operators q ± α are defined by their SU(2) commutation relations In order to model excitation relaxation time using the Hamiltonian H Q , we consider two separate approximation schemes, a fully integrable finite size model with 4 sites, and a dynamic Monte Carlo algorithm applied to an extended periodic boundary condition model with 200 sites.For the first approximation, using the short range nature of the Hamiltonian, we focus on a central vertex coupled to three isolated neighbors via the flip flop term of Eq. (S12).The result is a 4 4 dimensional Hilbert space system which we time evolve exactly.For the initial ensemble, the system is set at thermal equilibrium and projected into only states with a central q eigenvalue of q 0 = + 3 2 .We calculate the central excitation probability P (q 0 = + 3 2 ) as a function of time.The resulting time dependence of the q 0 = + 3 2 state shows a temperature independent relaxation time, in agreement with experimental results.The main qualitative change as a function of temperature is that the large time asymptote of the time average of P (q 0 = + 3 2 ) increases monotonically with temperature, as one would expect from ergodicity.These features can be seen in Fig. S8.
For the second approximate treatment of (S12), we utilize the dynamic Monte Carlo algorithm as presented in Ref. 4. The system we consider has dimensions of 10 unit cells by 10 unit cells, with 2 sites per unit cell, for a total of 200 sites, under periodic boundary conditions.We time evolve the system according to the dynamic Monte Carlo algorithm through thousands of hopping processes, recording any excitations that appear in the system and computing the average lifetime of stationary excitations.In this case, it is straightforward to model the system with the addition of a potential energy barrier to magnetic charge hopping.The experimental result of temperature independent relaxation time is reproduced for the case of this barrier being set to zero, as shown in Fig. 4b.Given the fact that exchange energy is expected to be much larger than the estimated barrier height, ∼ 20 K (as discussed below), the magnetic charge relaxation between neighboring vertices can be considered barrier-free.

B. Classical spin dynamics model
The large number of spins per honeycomb joint in the experiment, by conventional wisdom, would suggest the possibility of modeling the system according to the purely classical treatment of a Landau-Lifshitz-Gilbert (LLG) equation of motion.Here we provide such treatment while following the same philosophy as that above of restricting ourselves to a minimalistic model.As the starting point for the above kagome spin ice was the nearest neighbor spin exchange coupling, we start with the classical analogue of this term in our LLG Hamiltonian.We further posit an anisotropy term in order to favor spins pointing along the honeycomb joint axes described by the êi unit vectors.The result is the Hamiltonian where S i are normalized classical spin variables.Physically a term like K 1 arises in an effective field theory typically with dominant contributions from shape anisotropy; however this term can be used to wrap up a variety of both intrinsic and emergent anisotropic effects all into a single, phenomenological parameter.As such K 1 's optimal value can potentially be difficult to estimate both experimentally and theoretically.However, we have used the typical value for permalloy, K 1 ∼ −1 × 10 4 J/m 3 .In order to model the effects of temperature, we include a time dependent stochastic field B th in the LLG equations of motion, following the method of Ref. 5, which satisfies the time averaged expectation value relations B th (t) = 0, (S14) where γ is the gyromagnetic ratio, g is the dimensionless Gilbert damping parameter, M s is the saturation magnetization, V is the spatial volume of a magnetic unit cell, and δt is the time step size of the numerical integration.We numerically integrate a model consisting of three hexagons sharing three edges, with each edge corresponding to one spin S i (t), for a total of 15 spins.In order to model the dynamics of a classical version of a Q = 3 excitation, for the initial state, the three central spins are oriented outward, with the rest of the spins arranged in an energy minimizing orientation along their respective axes.For ferromagnetic coupling, this central Q = 3 excitation state corresponds to an unstable fixed point of the LLG equations of motions.For sufficiently small temperatures, corresponding to small thermal fluctuations B th (t), this model predicts an identical relaxation time for a given arrangement over a broad range of temperatures.However, we find that when the root mean square of B th (t) becomes comparable to the anisotropy field B an ∝ K 1 , the relaxation time begins to decrease monotonically as a function of temperature, as shown in Fig. 4a.
The existence of this anisotropy dependent crossover allows us to estimate a thermalization crossover temperature T c for the classical model.The values used for this estimate are δt = 10 −11 s, g = 0.2, M s = 6.6 × 10 5 A/m, and V = 264 nm 3 .The resulting crossover temperature estimate is on the order of T c ∼ 20 K, which is much too small to be consistent with the experimental results.The substantially better agreement of the experiment with quantum modeling provides evidence for our hypothesis on the robust quantum mechanical nature of the ultra-small artificial magnetic honeycomb lattice.
We again note that an effective value of K 1 can be difficult to estimate due to its emergent nature and due to the fact that the shape anisotropy intrinsically depends on the geometry of the material, which in this case is highly nontrivial.In spite of this, we believe that our theoretical model supports our conclusions of the dominant physics at play.

FIG. 2 :
FIG. 2:Neutron spin echo spectroscopy revealing dynamic properties in artificial honeycomb lattice without any external tuning parameter.a, Color map of NSE data on 6 nm thick honeycomb lattice with net neutron spin polarization (+Z) − (−Z).Bright localized scattering pixels are identified with q = 0.058 Å−1 and q = 0.029 Å−1 (only partially visible near x = 0 due to geometrical limit of the instrument).NSE data is obtained on a large parallel stack of 125 (117) honeycomb samples of 6 nm (8.5 nm) thickness.b,c, Spin echo profile at T = 300 K in 6 nm and 8.5 nm thick samples, respectively, at t = 0.02 ns neutron Fourier time.Strong signal-to-background ratio is detected in the NSE signal.d, The echo becomes weaker at higher Fourier time (t = 0.26 ns), indicating the dominance of charge defect qm relaxation at lower time scale.e,f, Spin echo profile at T = 4 K in 6 nm and 8.5 nm thick samples, respectively, at t = 0.02 ns neutron Fourier time.Self-propelled charge dynamics prevails at low temperature.

FIG. 3 :
FIG. 3: Temperature dependence of magnetic charge defect's relaxation time.a,b, S(q, t)/S(q, 0) versus neutron Fourier time at T = 300 K for the 6 nm sample (a) and the 8.5 nm sample (b).The magnetic charge relaxation time τm is obtained from an exponential fit of the data.c, τm versus T in the 6 nm sample.Charge defect qm relaxes at τ ∼ 20 ps timescale at T = 4 K without any external stimuli The plot also shows the temperature independence of τm.Inset: same as (a), but for T = 4 K.

FIG. 4 :
FIG. 4: Quantum mechanical nature of qm relaxation process.a,b, Temperature dependence of τm extracted from quantum (a) and classical (b) numerical simulations.Only the quantum Monte-Carlo simulations with negligible barrier height, compared to exchange constant, results in temperature independence of τm, as found in experimental data.It suggests effortless relaxation of qm, similar to quantum entity.
FIG.S1: Design of modified NSE instrument, utilized for the magnetic charge relaxation experiment.Inset shows the stack of samples sealed in an aluminum can, which serves as a π flipper.Additional magnetic coils are installed to enable the XYZ polarization analysis.

5 FIG
FIG. S1: Color detector image of the 6nm honeycomb lattice at room temperature with neutron's polarization Z+ and Z−.The color scale on the right shows the detected neutron intensity normalized by the proton charge.With the same color scale, while no meaningful spectral weight is detected in the Z− polarization, in the Z+ polarization case, clear spectral weight is observed in the detector region corresponding to q ∼ 0.058 Å−1 and q ∼ 0.029 Å−1 .The intense signal at the left edge of the detector is the direct beam as the instrument hits its geometrical constraint limit.Fig.2ais obtained by subtracting the Z− intensities from the Z+ intensities.

FIG. S2 :
FIG. S2: The q distribution and echo signals of the 6 nm honeycomb lattice measured at T = 4 K on the detector.Left panel shows the flux-weighted q distribution across the whole detector in ToF tbin 17, corresponding to a Fourier time of 0.02 ns, the numbers in the individual pixels denote the q value in the unit of 0.001 Å−1 .Right panel shows the echo signals in ToF tbin 17 on the detector.Every 2 pixels along X and Y directions are grouped for the clearness of the demonstration.The intense direct beam signals on the left edge of the detector are not shown.This image together with the color detector image in Fig. 2a guided the grouping of detector X, Y pixels in data analysis.

FIG. S2 :
FIG.S2: Color detector image of the 6nm honeycomb lattice at room temperature with neutron's polarization Z+ and Z−.The color scale on the right shows the detected neutron intensity normalized by the proton charge.With the same color scale, while no meaningful spectral weight is detected in the Z− polarization, in the Z+ polarization case, clear spectral weight is observed in the detector region corresponding to q ∼ 0.058 Å−1 and q ∼ 0.029 Å−1 .The intense signal at the left edge of the detector is the direct beam as the instrument hits its geometrical constraint limit.Fig.2ain the main text is obtained by subtracting the Z− intensities from the Z+ intensities.

5 FIG
FIG.S2: Color detector image of the 6nm honeycomb lattice at room temperature with neutron's polarization Z+ and Z−.The color scale on the right shows the detected neutron intensity normalized by the proton charge.With the same color scale, while no meaningful spectral weight is detected in the Z− polarization, in the Z+ polarization case, clear spectral weight is observed in the detector region corresponding to q ∼ 0.058 Å−1 and q ∼ 0.029 Å−1 .The intense signal at the left edge of the detector is the direct beam as the instrument hits its geometrical constraint limit.Fig.2ain the main text is obtained by subtracting the Z− intensities from the Z+ intensities.

FIG. S2 :
FIG. S2: The q distribution and echo signals of the 6 nm honeycomb lattice measured at T = 4 K on the detector.Left panel shows the flux-weighted q distribution across the whole detector in ToF tbin 17, corresponding to a Fourier time of 0.02 ns, the numbers in the individual pixels denote the q value in the unit of 0.001 Å−1 .Right panel shows the echo signals in ToF tbin 17 on the detector.Every 2 pixels along X and Y directions are grouped for the clearness of the demonstration.The intense direct beam signals on the left edge of the detector are not shown.This image together with the color detector image in Fig. 2a guided the grouping of detector X, Y pixels in data analysis.

FIG. S3 :
FIG. S3: The q distribution and echo signals of the 6 nm honeycomb lattice measured at T = 4 K on the detector.Left panel shows the flux-weighted q distribution across the whole detector in ToF tbin 17, corresponding to a Fourier time of 0.02 ns, the numbers in the individual pixels denote the q value in the unit of 0.001 Å−1 .Right panel shows the echo signals in ToF tbin 17 on the detector.Every 2 pixels along X and Y directions are grouped for the clearness of the demonstration.The intense direct beam signals on the left edge of the detector are not shown.This image together with the color detector image in Fig. 2a in the main text guided the grouping of detector X, Y pixels in data analysis.
FIG. S3: Echo profiles of the 6 nm honeycomb lattice at different Fourier times measured at room temperature, related to Fig. 2. The Fourier time of the echoes are labeled on the graphs.Clearly, as the Fourier time increases, the error bars become large and the echoes exhibit smaller amplitudes.

8 FIG
FIG.S5:The normalized intermediate scattering function of the 6 nm honeycomb lattice at q ∼ 0.058 Å−1 at

8 FIG
FIG.S5:The normalized intermediate scattering function of the 6 nm honeycomb lattice at q ∼ 0.058 Å−1 at

FIG. S6 :
FIG. S6: Measured paramagnetic resolution function at SNS-NSE and simulated elastic resolution function in soft matter regime.

8 FIG 7 FIG
FIG. S4: Echo profiles of the 6 nm honeycomb lattice at Fourier time 0.02 ns, measured at intermediate temperatures.Put together with Fig. 2b, and Fig. 3a-b, the persistence of a good signal-to-background ratio from room temperature to the lowest measurement temperature (4 K) demonstrates the high quality of the data.h=6nm FIG. S8: Exact time evolution of a four site quantum spin ice model ensemble with a central Q = 3 excitation initial state.The curves show the density of the central charge over time for three different initial temperatures of the surrounding sites.