Monopolar and dipolar relaxation in spin ice Ho2Ti2O7.

Ferromagnetically interacting Ising spins on the pyrochlore lattice of corner-sharing tetrahedra form a highly degenerate manifold of low-energy states. A spin flip relative to this "spin-ice" manifold can fractionalize into two oppositely charged magnetic monopoles with effective Coulomb interactions. To understand this process, we have probed the low-temperature magnetic response of spin ice to time-varying magnetic fields through stroboscopic neutron scattering and SQUID magnetometry on a new class of ultrapure Ho2Ti2O7 crystals. Covering almost 10 decades of time scales with atomic-scale spatial resolution, the experiments resolve apparent discrepancies between prior measurements on more disordered crystals and reveal a thermal crossover between distinct relaxation processes. Magnetic relaxation at low temperatures is associated with monopole motion through the spin-ice vacuum, while at elevated temperatures, relaxation occurs through reorientation of increasingly spin-like monopolar bound states. Spin fractionalization is thus directly manifest in the relaxation dynamics of spin ice.

INTRODUCTION

A rare-earth atom in a crystalline solid can form an isolated Ising-like doublet with large angular momentum ±ℏJ where low temperature reversal only occurs through transverse field–driven quantum tunneling. At elevated temperatures, the crystal field energy barrier separating reversed states can be surmounted over a time scale τ following the Arrhenius form (–)

The effective barrier Δ is near the first excited crystal field level, while τ0 [≃10−12 s for Ho2Ti2O7 ()] is the lifetime of the crystal field ground state. The divergence of τ(T) upon cooling complicates the attainment of thermal equilibrium among interacting Ising spins. When inter-spin interactions are unfrustrated and local order implies long-range order, adequate spin dynamics is typically maintained all the way to the critical temperature so that order can be achieved before spin tunneling ceases. Frustrated magnets, however, have reduced ordering temperatures and can support topologically protected defects that can become kinetically trapped and inhibit the formation of long-range order. The temperature dependence of the collective magnetic relaxation time reflects the development of short-range order and the nature and dynamics of topological defects that limit the correlation length. An illustrative example is found in quasi–one-dimensional SrHo2O4 where soliton defects in spin chains with frustrated interactions become trapped as spin order on the neighboring unfrustrated chains leads to an increase in the spin tunneling time on the frustrated chains ().

Here, we examine spin relaxation in classical spin ice (, ), where spins occupy the vertices of corner-sharing tetrahedra (Fig. 1A). Ferromagnetic interactions (Jeff) define a degenerate manifold of short-range ordered states satisfying the magnetic analog of the Bernal-Fowler ice rule () wherein as many spins point in as point out of each tetrahedron (Fig. 1B) (, ). While theory predicts long-range order at sufficiently low T (), AC susceptibility measurements in Dy2Ti2O7 show that spin relaxation times increase beyond the single-spin quantum tunneling time scale of ∼10−4 s upon cooling below 2 K and exceed 105 s for T < 0.3 K (–). As a result, the spin system falls out of thermal equilibrium for any realizable cooling protocol and becomes trapped in a spin ice state with near-perfect local order but no long-range order. Explorations of spin dynamics in Ho2Ti2O7 are less complete but reveal qualitatively similar results (–).

Fig. 1

Spin ice state and time-resolved diffuse neutron scattering from Ho2Ti2O7.

(A) Cubic pyrochlore lattice with Ho (blue) on corner-sharing tetrahedra. Our (HHL) scattering plane is indicated. (B to D) Spin configurations with and without monopoles. The spin-ice manifold has two spins pointing in and two pointing out of every tetrahedron (B). A spin flip (B to C) creates a pair of “three-in-one-out” (red cloud) and “three-out-one-in” (blue cloud) defects. No further ice rule violation is needed to separate monopoles (C to D). A monopole-antimonopole pair can form a dipole (C). (E to I) Elastic neutron scattering in the (HHL) plane at T = 0.95 K while being subjected to 10-s, 250-Oe field pulses [inset in (F)] applied along [arrow in (D)]. While zero-field diffuse scattering (E) indicates the spin-ice ground state, the field-induced difference (F) concentrates at the Bragg points within white frames. The time-dependent intensity near (002) is unfolded along the (HH0) (G) and (00L) (H) directions. The integrated intensity ∫ ΔI(Q, t)d2Q (I) documents the opposite time dependence of the net Bragg intensity at (002) and (220) versus the net non-Bragg intensity (blue). The red line is the best fit using Eq. 7 (Materials and Methods), while the blue line is its mirror image relative to the x axis. For (E to I), zero-field and 250-Oe conditions refer to the 12s < Δt < 20s and 2s < Δt < 10s time windows, respectively [inset in (F)]. Error bars represent ±σ counting statistics.

The topological defects in spin ice are magnetic monopoles that result from spin fractionalization (Fig. 1, C and D) (, ) and are experimentally supported by neutron scattering experiments (, ). Monte Carlo simulation and Debye-Hückel theory based on a Coulomb gas model have achieved great success in accounting for heat capacity (, ) and AC susceptibility measurements in the higher temperature regime (, –). However, the sharp rise of the spin relaxation time deep in the spin ice phase cannot be captured without using a temperature-dependent metropolis time (, , , ). At lower temperatures, the question as to whether monopoles become free because of entropic repulsion or remain as bound pairs because of Coulomb attraction () is not firmly resolved by existing AC susceptibility measurements.

Here, we report the discovery of a distinct thermal crossover in the magnetic relaxation dynamics of Ho2Ti2O7 spin ice from a low temperature regime with a well-defined Debye-like relaxation to a higher temperature regime with a broader and faster relaxation spectrum. Enabled by a new class of ultrapure Ho2Ti2O7 single crystals (), a novel time-resolved magnetic neutron scattering technique, and a broadband AC magnetometry, our experiments resolve apparent discrepancies in the literature where experiments covering different temperature ranges in more disordered samples probed one or the other but never both of these distinct relaxation regimes (, ). The totality of the data allows association of the two regimes with the distinct dynamics of monopoles and dipoles, respectively.

RESULTS

Time-dependent neutron scattering

Figure 1E shows the momentum (Q)–dependent elastic magnetic neutron scattering S(Q) in the (HHL) scattering plane deep in the spin-ice regime of Ho2Ti2O7 (T = 0.95 K). Devoid of Bragg peaks, the exclusively diffuse magnetic scattering surrounding high-symmetry points on the Brillouin zone boundary, such as (001), (003), and , reflects local noncollinear ferromagnetic Ising spin order on each tetrahedron (Fig. 1B). The faint “bow tie”–like pinch points at the (002) and (111) zone centers indicate the divergence-free nature of magnetization in the spin-ice manifold (, ).

To probe the magnetic relaxation dynamics of spin ice, we measured the time dependence of magnetic neutron diffraction from spin ice while being exposed to a square pulse sequence of magnetic field alternating between 250 Oe and zero field with a 20-s period (Fig. 1F, inset). By registering the time of arrival of elastically scattered neutrons relative to the field pulse, we obtained the time and momentum dependence of the neutron scattering cross section. This novel method is able to probe magnetic dynamics on time scales down to 0.1 ms and is detailed in Materials and Methods.

Figure 1F shows the impact on scattering of applying H = 250 Oe along and perpendicular to the (HHL) scattering plane. Resolution-limited magnetic Bragg peaks are clearly observed at (002) and (220), where nuclear Bragg diffraction is forbidden or weak, with magnetic diffraction intensity proportional to the induced magnetization squared. There is also a general reduction in diffuse scattering throughout the (HHL) plane that is Q independent, uncorrelated with S(Q) to within statistical accuracy, and not associated with changes in the intensity of pinch point scattering (see note S1 and fig. S1). These observations negate field-induced changes in the monopole density and imply that the net magnetization is induced through motion of magnetic monopoles along the field direction (Fig. 1, C and D). Changes in the Q dependence of S(Q) for Q⊥H were previously shown to require considerably larger magnetic fields ().

Figure 1G to I show the time dependence of the coherent and of the diffuse magnetic scattering driven by the time-dependent applied field shown in Fig. 1F. The resolution-limited nature of the magnetic Bragg peak indicates that the induced magnetization is homogeneous throughout the coherence volume of the scattering experiment, which extends over (~500 Å)2 perpendicular to the applied field but just 9 Å along the field direction (because of the vertically focusing monochromator and analyzer). The integrated magnetic Bragg intensity at Q = (002) (Fig. 1I) is a measure of the time-dependent transverse-to-Q magnetization squared . For comparison and in the same units, we present the time dependence of the diffuse magnetic scattering throughout the accessible region in Fig. 1F excluding Bragg peaks. As anticipated from the total moment sum rule, what is gained in magnetic Bragg scattering while the field pulse is on is lost in diffuse scattering that is independent of Q⊥H. The contrast between the Q independence of the time-dependent diffuse scattering and the sharp Q dependence of the time-dependent Bragg scattering that develops from it indicates that the magnetization arises from uncorrelated point-like objects whence monopoles that remain uncorrelated within the plane perpendicular to the applied field.

Having established that time-dependent Bragg intensity at (002) is a measure of and, hence, monopole displacement along H ⊥ (002), we can extract the time scale over which monopole drift ceases following a step change in magnetic field by analyzing I(t) datasets such as those shown in Figs. 1I and 2A. Upon cooling from 1.3 K to 0.6 K, this time scale increases from milliseconds to hours. A stretched exponential description of M(t) (see Materials and Methods) provides a good account of the intensity data (solid lines in Figs. 1I and 2A), and the corresponding T-dependent time constant τ and exponent β are shown as solid black symbols in Fig. 3B and 3C, respectively.

Fig. 2

Time-dependent neutron diffraction and frequency-dependent susceptibility, spanning more than nine decades of time scales.

(A) Time-dependent diffraction intensity at (002) under periodic perturbation by H = 100-Oe field pulses for T= 0.6 K to 1.3 K, scaled to the expected intensity under equilibrium field-on conditions (Materials and Methods). These time sequences represent responses to step changes in the magnetic field (Fig. 1F, inset), with the pump period varying from 10−2 to 104 s. Solid lines represent best fits using Eq. 7. (B and C) AC susceptibility χ(f) measured by a superconducting quantum interference device (SQUID) magnetometer in a dilution fridge down to 1 mHz (B) (Materials and Methods) and by a Quantum Design (QD) ACDR (AC susceptometry in a dilution refrigerator) down to 10 Hz (C) (Materials and Methods). Real (χ′) and imaginary (χ″) parts of χ(f) are in the top and bottom panels, respectively. Corrections for the response function, demagnetization effects, and calibration to absolute units were applied as described in Materials and Methods. Only selected curves from our ACDR measurements are displayed (all curves are presented in fig. S2). Solid lines in (B) represent fits using two modes (Eq. 3) with β = 0.728 and β = 0.492. In (C), solid lines represent best fits to the single-mode Havriliak-Negami form (Materials and Methods) for the best guide to the eye. The discrepancy at low frequency indicates the second mode that is resolved in the lower-frequency SQUID data (B). Across all panels, colors consistently indicate temperature (see color bar). Magnetic fields are applied along the direction (arrow in Fig. 1D). The magnetic susceptibility is presented in electromagnetic units (emu).

Fig. 3

Two relaxation modes in Ho2Ti2O7.

(A and B) Temperature dependence of χ0 and τ for χ (blue) and χ (red) in Eq. 3 with β = 0.728 and β = 0.492. Filled circles and triangles are SQUID (Fig. 2B) and ACDR (Fig. 2C) measurements, respectively. Blue and red dashed lines in (A) guide the eye, while the black dashed-dotted line is χ0 = χ0 + χ0. (B) Time constants from neutron scattering (Fig. 2A), from our susceptibility measurements (Fig. 2, B and C), and from (, ). Solid lines are best fits to the Arrhenius law with parameters in Table 1. τ from neutron scattering (black points) is three times longer than the averaged τ and τ. This may result from imperfect demagnetization correction (Materials and Methods). A plot of logτ versus 1/T is in the inset. (C and D) Line shape analysis of neutron and susceptibility data below T ~ 1.3 K using a single-mode form. The T dependence of β is apparent from susceptibility (Eq. 2) and neutron scattering (Eq. 7) and the half-width half-magnitude (HWHM) of χ″[log (f)] on the low-frequency (HWHM−; square) and high-frequency (HWHM+; circle) sides. The gray dashed line in (D) indicates the expected HWHM of 0.57 log(Hz) for Debye relaxation. The T-independent asymmetry indicated by HWHM+ > HWHM− below 0.9 K in our data (D) is consistent with β = 0.728(4) < 1 in (C). The increase in HWHM+ for T > 0.9 K indicates the onset of the faster mode. The peak widths from () [(D), open symbols] demonstrate ∼50% broadening and a lower onset temperature for the high-T mode in the more disordered floating zone (FZ) samples.

Table 1

Arrhenius parameters for Ho2Ti2O7 and Dy2Ti2O7.

The Arrhenius form τ(T) = τ0 exp (Δ/T) (Eq. 1) is used to fit the temperature dependence of τ and τ in Ho2Ti2O7. The frequency and temperature range of the instruments used indicates that τ is probed in () and τ is probed in (). The samples investigated in this work are from the TSFZ method (Materials and Methods) (), while the samples studied in (, ) are from traditional FZ growth (Materials and Methods). Arrhenius parameters for Dy2Ti2O7 () and Jeff values for both Ho2Ti2O7 and Dy2Ti2O7 (, ) are provided for comparison. Arrhenius fitting of τ(T) from neutron measurements gives Δ= 18(2) K and τ0 = 4(8) × 10−8 s, without distinguishing τ and τ (Materials and Methods).

Ho2Ti2O7, Jeff ≃1.8 K
	τ0l (s)	τ0h (s)	Δl (K)	Δl/Jeff	Δh (K)	Δh/Jeff
Currentstudy	2.3(8) ×10−7	4(1) ×10−9	15.3(3)	8.5	18.3(3)	10.2
(14, 15)	4.5(9) × 10−7	9.5(7) × 10−11	10.4(2)	5.8	17.5(9)	9.7
Dy2Ti2O7, Jeff ≃ 1.1 K
	τ0 (s)		Δ (K)		Δ/Jeff
(12)	4.06 × 10−7		9.79		8.9

Frequency-dependent magnetic susceptibility

Sampling the magnetic relaxation in the frequency domain provides important complementary information, so we carried out AC susceptibility measurements, which also extend to higher frequencies and temperatures than the neutron scattering data (Fig. 2, B and C, and fig. S2). Unexpectedly, two distinct characteristic frequencies are apparent for T = 1.5 K (Fig. 2B). This contrasts with the well-defined single-mode response that we observe at low temperatures and that is reported in the literature (, ). As χ(f) at low temperature can be described by the empirical Cole-Davidson form (see Materials and Methods and fig. S3)we use a superposition of two such response functions to describe the bimodal spectrum at T = 1.5 K (Fig. 2B)

Here, χ and χ refer to the response that dominates at low and high temperatures, respectively. χ0 and χ0 are the corresponding static susceptibilities, and τ and τ are the characteristic relaxation times. β and β characterize the distribution of relaxation time scales. As for the stretched exponential in the time domain, β = 1 represents Debye relaxation with a single time scale, while β < 1 describes an asymmetric spectral function with extra weight in a high-frequency tail (). To stabilize the fitting analysis, we fix β and β to be temperature independent and determined them in a simultaneous global fit to all the susceptibility data, which yields β = 0.728(4) and β = 0.492(8) (for detailed analysis, see Materials and Methods and fig. S4). This suggests distinct microscopic characters for the l and h components of the relaxation response.

Figure 3B also includes, as open symbols, data from two separate previous studies of the AC susceptibility of Ho2Ti2O7 (, ) that initially appear to be inconsistent with each other and with the present results. We attribute the three to four decades longer time scales in the current work (Fig. 3B) to our highly stoichiometric traveling solvent floating zone (TSFZ)–grown crystals (see also Materials and Methods) (). Oxygen vacancies and Ho3+ stuffing associated with conventional floating zone (FZ)–grown samples may offer nucleation centers for monopole creation () and effective transverse fields on the non-Kramers Ho3+ ions that increase the relaxation rate ().

Two modes were not previously resolved in Ho2Ti2O7. This may, in part, be because the frequency where the two modes can be clearly separated is reduced by two orders of magnitude and becomes more experimentally accessible in higher-quality samples (Fig. 3B). Nonetheless, there are, in fact, indications of a faster relaxation process in the previously published line shape analysis for χ(f) in lower-quality FZ samples. Specifically, () documents the abrupt appearance of an asymmetric line shape upon heating, which is also apparent in a single-mode analysis of our data for high-quality TSFZ samples (Fig. 3, C and D).

Further distinguishing the two modes, the inset to Fig. 3B shows that disorder affects the high- and low-T relaxation processes differently. For the high-T process (red symbols), the Arrhenius lines have identical slopes, indicating that impurities simply reduce the asymptotic relaxation time τ0 (Table 1). For the low-T process (blue symbols), it is, however, the slope of the Arrhenius line and, hence, the activation energy Δ that is reduced from 15.3(3) K to 10.4(2) K by impurities (Table 1).

DISCUSSION

Monopole dynamics at low temperatures

The Q⊥-independent nature of the field-dependent diffuse magnetic scattering (Fig. 1, F and I; note S1; and fig. S1) indicates that magnetic relaxation in the low-T regime for T < 1.25 K is driven by the motion of thermally activated monopoles along the applied field direction. In spin ice, magnetic monopoles can be mapped to magnetic charges that are connected through a chain of spin dipoles known as “Dirac strings” (, , , ). Within this dumbbell model, the monopole conductivity can be written as , where μ0 is the vacuum permeability and μ(T), Qm, and n(T) are the monopole mobility, charge, and number density, respectively. The response of magnetic monopoles to an external magnetic field thus should follow the magnetic version of “Ohm’s law”: Jm = σmH (), which implies a magnetization relaxation time . Charge relaxation in a semiconductor where τe = ϵ0/σ provides the electrical analog. In this framework, the rise in the relaxation time upon cooling (Fig. 3B) simply indicates the increasing monopole resistivity anticipated for a monopole “semiconductor.” Much as the resistivity of a semiconductor, the increase in the magnetization relaxation times τ(T) and τ(T) upon cooling can be described by the Arrhenius law (Eq. 1).

The inferred barrier height Δ = 15.8 K (Table 1) exceeds the energy of a free monopole Δmono, which is estimated to be 5.7 K for Ho2Ti2O7 and 4.35 K for Dy2Ti2O7 through Debye-Hückel theory and Monte Carlo simulation of the Coulomb gas (, ). This indicates that the monopole mobility μ(T) decreases with T in proportion to a power of the monopole density: μ ∝ nη with η ≃ 1.8 for high-quality Ho2Ti2O7 and η ≃ 1.2 for Dy2Ti2O7 (Table 1) (, ). Monopole motion is ultimately associated with quantum tunneling of flippable spins adjacent to monopoles with a tunneling rate determined by transverse fields from other spins or monopoles (–). Although it remains unclear how η can be directly associated with microscopic parameters of the spin ice model, different values of η for Ho2Ti2O7 and Dy2Ti2O7 are not surprising, given the distinct non-Kramers versus Kramers nature of Ho3+ versus Dy3+ and the different relative strength of exchange interactions and dipolar interactions in the two compounds (, ). Hyperfine coupling to nuclear spins may also play a role and is different for Dy3+ and Ho3+ (). None of these effects are captured by the Monte Carlo simulation because they affect the metropolis time. Neglecting the T dependence of the metropolis time, the Monte Carlo simulation yields μ ∝ 1/T (), which is inconsistent with the data as previously noted (, , , ).

The width of the peak in χ″[ log (f)] is a measure of the distribution of relaxation times. In glass-forming systems, the relaxation time distribution typically broadens upon cooling (). A concomitant decrease in β is an alternate indicator of glassy heterogeneous relaxation as in proton glasses, which display logarithmic time dependence (β → 0) (). In Ho2Ti2O7, on the other hand, the width of the peak in χ″[ log (f)] actually narrows upon cooling. Figure 3D shows that the half-width half-magnitude (HWHM) of the peak in χ″[ log (f)] approaches the limit associated with a single characteristic relaxation time scale. Correspondingly, our stroboscopic neutron scattering data show that β remains close to 1 down to T = 0.6 K (Fig. 3C) while τ increases by six orders of magnitude (Fig. 3B). This observation of Debye relaxation over a 105-s time scale in the low-T limit is unique to our knowledge. It indicates a single dominant relaxation time scale and process as might be associated with the independent motion of monopoles at low density through the homogeneous spin-ice vacuum. Independent monopole motion in the low-T regime is consistent with Ho2Ti2O7 not supporting monopole crystallization (, ). This interpretation of the low-T magnetic relaxation process is supported by the deviation of τ(T) from the Arrhenius form upon heating (Fig. 3B), which indicates loss of homogeneity as the “ice rule” is violated at elevated monopole densities.

Relaxation through spin dipoles at high temperatures

We now turn to the T > 1.2 K regime where a distinct high-T relaxation process emerges (Fig. 3A). Monopoles subject to competing hopping processes with time scales given by τ(T) and τ(T), respectively, would result in a branching ratio χ0/χ0 = (τ0/τ0) exp [(Δ − Δ)/T], implying a crossover where χ0 ≈ χ0 for T = (Δ − Δ)/ ln (τ0/τ0)= 0.69 K. This is inconsistent with the observed crossover temperature of T ≃ 1.20 K (Fig. 3A). If we instead associate the high-T mode with an activation energy scale ΔE, then we obtain χ0/χ0 = C exp (ΔE/T), which fits the data in Fig. 3A with a numerical constant C = 3.4(6) × 10−4 and ΔE = 9.7(9) K (fig. S5A). Both the high crossover temperature and the small numerical factor C < < 1 indicate that the high-T mode arises from a magnetizable entity that is distinct from isolated monopoles. This conclusion is supported by the very different values of β and τ0 for the two modes (Table 1). The activation energy ΔE= 9.7 K is close to the energy cost to flip a spin out of the spin-ice manifold (, ), so the Bjerrum monopole-antimonopole pair is a natural candidate (, ). As its characteristic dimension decreases on warming, the Bjerrum pair becomes a spin flip relative (Fig. 1C) to the spin-ice manifold when lB(T) approaches the diamond lattice constant ad. The corresponding magnetic susceptibility is , where nB is the pair density (note S2). Associating χ0(T) with χB(T), we can infer nB(T), which is qualitatively consistent with the prediction from Debye-Hückel theory (, ), although the experiment yields a sharper increase in nB(T) with T and saturation is achieved at lower temperature (fig. S5B). The lower value of β = 0.492(8) (fig. S4) indicates a broad spectrum of relaxation times, which may arise from inhomogeneity and short-range interactions as in glass-forming systems (, , ). In addition, the asymptotic microscopic tunneling time scale τ0 (Eq. 1) decreases by a factor of 20 in the more disordered crystals (Table 1). Both features are consistent with the greater complexity of Bjerrum pair reorientation as compared to individual monopole motion.

Coexistence of magnetic monopoles and Bjerrum bound pairs was previously reported through real-time imaging of artificial spin ice (). While monopolar and dipolar entities may not be experimentally distinguishable in heat capacity and static magnetization data, the dynamic crossover reported here indicates that they may be distinguished in magnetic relaxation dynamics. While a crossover in relaxation dynamics is not apparent in published susceptibility data for Dy2Ti2O7 (fig. S6) (, ), Dy2Ti2O7 differs from Ho2Ti2O7 in the nature of the single-ion spin and their interactions and in the effects of disorder. The approach to slow Debye relaxation in the low-T regime of Ho2Ti2O7 is inconsistent with conventional spin freezing but is evidence of relaxation through motion of a diminishing density of monopoles through the spin-ice manifold. That this regime is realized in high-quality single crystals is encouraging for the prospects of coherent quantum dynamics of monopoles () in quantum siblings such as Ce2Zr2O7 (, ).

MATERIALS AND METHODS

Samples and demagnetization factors

The growth technique and detailed characterization of the Ho2Ti2O7 samples investigated in this work can be found in () and are briefly summarized here. All the samples were cut from the same single-crystal rod grown by the TSFZ technique. The crystals investigated in (, ) were grown at higher temperatures using the FZ method. As reported in (), refinement of x-ray diffraction data revealed off-stoichiometry corresponding to x > 0.01 in Ho2-(Ti2Ho)O7-δ for the FZ-grown crystals, which also results in a larger lattice parameter. AC susceptibility data show that the FZ samples have notably higher relaxation rates than TSFZ samples. TSFZ crystals are also more transparent than FZ crystals. The dimensions of each sample and the demagnetization factors (, ) used in the data analysis are listed in Table 2.

Table 2

List of Ho2Ti2O7 crystals used in the experiments.

The table provides the dimensions of each sample used for neutron scattering, SQUID magnetometry, ACDR, and with the AC magnetization measurement system (ACMS) (Materials and Methods). The corresponding demagnetization factors listed in the third column were calculated following (, ).

Experimentalmethods	Dimension (mm) (l isalong the fielddirection)	Demagnetizationfactor
Neutron scattering	Cylinderdiameter = 7.0, l = 8.0	0.28
SQUID	Cross section = 3.0 ×0.53, l = 5.0	0.09
ACDR and ACMS	Cross section = 0.92 ×0.64, l = 2.0	0.16

Time-resolved neutron scattering

The neutron scattering experiment was performed in a pumped Helium-3 cryostat. At each temperature, a magnetic field was periodically turned on and off, with dwell time in each state set to be always longer than, or comparable to, the relaxation time scale. Neutrons scattered from the sample were recorded in a time series referenced to the most recent field on time. In approaching a given temperature for measurements, no systematic control of cooling rates was imposed. We verified at several temperature points that data acquired after different thermal cycling histories are consistent. Note that certain commercial equipment, instruments, materials, or software are identified in this section to foster understanding. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the materials or equipment identified are necessarily the best available for the purpose.

With GE varnish at the interface, the Ho2Ti2O7 crystal was fastened to the top of a single-crystal sapphire rod using Teflon strips. The rod was mounted in a holder made of oxygen-free high–thermal conductivity copper affixed by GE varnish and black epoxy (STYCAST 2850FT). The copper holder was then threaded to the Helium-3 pot of the cryostat. Two RuO2 thermometers (Lakeshore Rox-102A) were separately mounted on the sapphire rod and to the top of the Ho2Ti2O7 crystal and read by the two channels of a Lakeshore 340 temperature controller. The magnetic field was generated by an aluminum solenoid, centered around the sample, and mounted around the inner vacuum chamber, which was made from single-crystalline silicon. A current of 1 A to 2 A was provided by a Dynatronix-DCR DPR20-15-30(XR) power supply with programmable output time sequences. Pump periods from 80 ms to 200 s were controlled by the internal clock of the power supply. Longer pump periods were externally controlled by a pulse generator (Berkeley Nucleonics Corporation Model 565).

Elastic neutron scattering was performed with incident energy Ei = 5 meV on the triple-axis spectrometers SPINS (Spin Polarized Inelastic Neutron Spectrometer) and MACS (The Multi Axis Crystal Spectrometer) at the NCNR (National Institute of Standards and Technology Center for Neutron Research) for (002) diffraction and diffuse scattering, respectively. Time series of events from the neutron detector(s) and the step changes of the magnetic field were recorded by a time stamper (General Electric, RS-DCS-107M4968) and then folded into a single pump period with t = 0 defined as the field on time. While measuring diffuse scattering at MACS, the sample rotation angle A3 was continuously swept through a 180° range, leading to a time-dependent coverage of Q space, asynchronous with the sequence of field steps. A coverage function in the two-dimensional A3-Δt space was accumulated to convert the event stream into scattering intensity with error bars in units of counts per second (Fig. 1, E to I, and fig. S1). Data from MACS were analyzed using Mslice in DAVE () with the event mode package updated on 4 May 2020. The time resolution of this neutron scattering method is determined by the distribution of time taken for a neutron to be detected following an elastic scattering event at the sample, i.e., , with l as the distance from sample to the position of detection, v as the neutron speed, and dtdet ≃ 3 μs as the pulse width for a neutron detection event. In our experiment, l ~ 1.2 m and , where the transverse and the longitudinal arise from the 1° beam divergence for diffraction and the diameter (depth) of the detector d = 25.4 mm, respectively. The speed of neutrons with incident energy of Ei = 5 ± 0.2 meV is v = 977 m/s with ~ 0.02. Therefore, dt ≈ 35 μs in our experimental configuration. On the pump side, the step change of current output from the power supply had a width of 50 μs and 150 μs for the rising and falling edges, respectively, which set an additional limit for time resolution. All combined, we estimate the time resolution of our pump-probe setup to be 10−4 s.

Frequency-dependent susceptibility measurements

AC susceptibility in absolute unit [electromagnetic unit (emu) Oe−1 cm−3] was measured by the commercial AC magnetization measurement system (ACMS) of the Quantum Design (QD) Physical Property Measurement (PPMS) for temperatures between 1.9 K and 7 K and frequencies from 10 to 104 Hz (fig. S2A). These data were used to normalize the lower-temperature measurements from 1 K to 3.7 K on the same sample using the AC susceptometry in a dilution refrigerator (ACDR) option of QD PPMS (Fig. 2C and fig. S2B) over the same frequency range. A superconducting quantum interference device (SQUID)–based magnetometer (iMAG 303 Multi-Channel SQUID System) with a 5-kHz low-pass filter was used to measure AC susceptibility down to f = 1 mHz in a dilution fridge (Fig. 2B). The sample was glued with GE varnish onto a silver plate that was pressed and held by a screw into the copper mount for the coil set. This mount, carrying the drive and pick-up coils, was anchored on the mixing chamber. A lock-in amplifier (Stanford Research, SR830) was used to provide a sinusoidal drive signal and pick up the in-phase and out-of-phase signals from the output of the SQUID control unit, as measures of χ′ (f) and χ″ (f), respectively. The remnant uncompensated signal, generated by the mismatched inductance of coils mounted on the mixing chamber, produced a temperature-independent background, which was measured at T = 0.5 K where the response from Ho2Ti2O7 is known to be too slow to generate any signal in the measured frequency range. The response function of the SQUID circuits was determined by comparing frequency scans at T = 4 K that were separately measured by SQUID and ACMS, which, meanwhile, calibrated the conversion factor from millivolt to emu for the SQUID magnetometer to be 3.1 × 10−7 emu/mV. All the presented AC susceptibility curves (Fig. 2, B and C, and figs. S2 and S3) were corrected for demagnetization effects using the equations given in () for centimeter-gram-second (cgs) system of units. Measurements from the SQUID and the ACDR match within 10%, as demonstrated by the χ(f) curves at T = 1.5 K (fig. S2C), after their independent calibration by ACMS measurements in the separate ranges of overlap.

Data analysis in time and frequency domains

All the least-squares fitting analysis was performed using the Global Optimization Toolbox of MATLAB, looking for convergence across 60 different starting points in parameter space. AC susceptibility measured in the frequency domain can be empirically described by the Havriliak-Negami form of (, , ), which reduces to the Cole-Davidson form (Eq. 2) with α = 1 and further to the Debye form when α = β = 1. We verified that the Cole-Davidson form is the minimal model to describe the measured χ(f) curves (fig. S3). This is consistent with previous investigations of Dy2Ti2O7, which argued that β ~ 0.7 to 0.8 is the ideal case of β = 1 modified by boundary effects (). The forms given by Eqs. 2 and 3 were used for single- and two-mode analysis, respectively. The simultaneous fitting of the real and imaginary parts of χ(f) automatically respects the Kramers-Kronig relationship.

We first describe the single-mode analysis whose ultimate failure at high temperature is clearly demonstrated in fig. S3. The Fourier transform of χ(f) given in Eq. 2 was used to describe the time-domain relaxation function M(t), in response to a step change of magnetic fieldwhere M0 is the equilibrium magnetization and is the three-parameter Mittag-Leffler function (, ), which was numerically evaluated (). For the demagnetization correction, we used the numerically generated M(t) array fromwhererepresents the response to a δ function impulse, which is a direct Fourier transform of χ(f) given by Eq. 2 (, ). At a given temperature, χ0 was extracted from susceptibility measurements, H0 is the external magnetic field, and D is the demagnetization factor.

The neutron scattering intensity at Q = 0 positions is related to M(t) throughwhere I0 is the nuclear scattering background and A is a numerical factor determined by the scattering geometry and cross section. When equilibrium is reached for t > > τ, the fully developed magnetic moment leads to neutron scattering intensity . The time-dependent curves presented in Fig. 2A are scaled as (I(t) − I0)/(Ieq − I0).

We note that Eq. 4 using the Mittag-Leffler function is analogous to the stretched exponential function (, ) but with the advantage that τ and β have the same physical meaning as in Eq. 2, which facilitates comparison to frequency domain measurements (Fig. 3, B and C). Despite the much larger field used for neutron scattering compared with the AC susceptibility measurements, we found that the field dependence of τ measured by neutron scattering is within 15% from 50 Oe to 150 Oe and can be understood as a result of pulsed field-induced heating.

For Ho2Ti2O7, use of this single-mode formalism leads to temperature dependence of the inferred value of β for temperatures above 0.9 K and discrepancy between neutron and susceptibility results at high temperatures (Fig. 3C). This can be understood as a consequence of the modified line shape of χ(f) resulting from the high-T mode, which affects β refined from neutron and susceptibility measurements differently because of their linear and logarithmic sampling in the time and frequency domain, respectively. For the same reason, applying two-mode analysis to time-dependent neutron scattering is not as accurate. Thus, we report the single time scale τ (black symbols in Fig. 3B) that mixes the contributions from τ and τ. The discrepancy between τ and the average of τ and τ measured by AC susceptibility can be attributed to the imperfect demagnetization correction. Demagnetization effects are complicated for the neutron scattering experiment because of the large sample and the large magnetic moment of Ho3+.

When applying the two-mode model (Eq. 3) to χ(f) measured by SQUID magnetometry (Fig. 2B), fitting at temperatures above T = 1.14 K converged without any constraints or particular choice of starting points. However, when the two characteristic frequencies get closer and the spectral weight of the high-T mode gets weaker at lower temperatures, the fitting analysis is ambiguous. In particular, the spectral weight at high frequency can be accounted for by β < 1. This is a direct consequence of the Cole-Davidson functional form (Eq. 2), which complicates the determination of χ0, β, and τ. In the following, we detail and legitimate the constraints used in our fitting analysis.

We noticed that β of the low-T mode is essentially the temperature-independent β below T ~ 0.9 K in the single-mode description (Fig. 3C), while β was found to fluctuate within 0.48 ± 0.03 with no discernible temperature dependence for T > 1.14 K (fig. S4A) when two modes can be unambiguously determined. Therefore, we applied a global optimization to reduce the redundant freedom and obtained β = 0.728(4) and β = 0.492(8) for all temperatures (fig. S4, B and C). We further found that convergence is facilitated by using starting parameters for τ in accordance with the Arrhenius form. These additional conditions were used for temperatures below 1.14 K. Fitting results for χ0 and τ with the above conditions are presented in Fig. 3 (A and B, respectively). The emergence of the high-T mode around T ≃ 0.9 K captured by χ0 (Fig. 3A) is consistent with signatures from the line shape (Fig. 3, C and D).

Measurements from the ACDR were fit solely by the high-T mode χ since the spectral weight of the low-T mode lies beyond the frequency window of the instrument (10 Hz to 104 Hz). The limited frequency range affects the overall fitting quality because the tails of the low-T mode create a background shape that escapes accurate modeling (Fig. 2C and fig. S2).

Least-squares fitting analysis

We denote x as the Nd independent variables associated with an experimental observable y(x) with standard deviation (SD) σ(x) described by an objective function f(x; p). p are the Np parameters that are constrained by the data to within an SD σ. The χ2 is given by

The fit result p = p is determined through minimizing χ2, i.e., , and the 1σ contour in the parameter p space was extracted through ND = Nd − Np is the number of degrees of freedom.