Sleuthing out exotic quantum spin liquidity in the pyrochlore magnet Ce2Zr2O7

The search for quantum spin liquids—topological magnets with fractionalized excitations—has been a central theme in condensed matter and materials physics. Despite numerous theoretical proposals, connecting experiment with detailed theory exhibiting a robust quantum spin liquid has remained a central challenge. Here, focusing on the strongly spin-orbit coupled effective S = 1/2 pyrochlore magnet Ce2Zr2O7, we analyze recent thermodynamic and neutron-scattering experiments, to identify a microscopic effective Hamiltonian through a combination of finite temperature Lanczos, Monte Carlo, and analytical spin dynamics calculations. Its parameter values suggest the existence of an exotic phase, a π-flux U(1) quantum spin liquid. Intriguingly, the octupolar nature of the moments makes them less prone to be affected by magnetic disorder, while also hiding some otherwise characteristic signatures from neutrons, making this spin liquid arguably more stable than its more conventional counterparts.


INTRODUCTION
The search for quantum spin liquids (QSL)-exotic forms of matter with fractionalized excitations-has drawn much attention of the physics community. While theories are no longer in short supply, tracking down materials has turned out to be remarkably tricky, in large part because of the difficulty to diagnose experimentally a state with only topological, rather than conventional, forms of order. Pyrochlore systems have proven a particularly promising ground, with the spin ices Dy/Ho 2 Ti 2 O 7 1 hosting a classical Coulomb phase [2][3][4][5] , with subsequent proposals of candidate QSLs in other pyrochlores 6 . Quantum fluctuations have a tendency to pick an ordered state from the classically degenerate ice manifold via the "order by disorder" 7 mechanism, making QSLs extremely rare in nature. However, it is known that certain QSLs exist as a robust phase in realistic theoretical models on the pyrochlore lattice 8 ; and one thus hopes that a material can be found that realizes a "quantum spin ice".
The recently discovered pyrochlore Ce 2 Zr 2 O 7 9,10 shows tremendous promise as the long sought-after QSL/quantum spin ice. Its neutron-scattering response shows a complete lack of magnetic ordering and the presence of a continuum of spectral weight 9,10 , a characteristic hallmark of fractionalized excitations. Additionally, no signature of an ordering transition at low temperature has been observed in the specific heat. These observations immediately led to the proposals that Ce 2 Zr 2 O 7 might provide a realization of a QSL.
What makes the cerium-based Ce 2 Zr 2 O 7 pyrochlore distinct from the earlier studied Yb-based quantum spin-ice candidate Yb 2 Ti 2 O 7 11-14 is the dipolar-octupolar nature [15][16][17] of the ground state doublet of the Ce 3+ ion, shown schematically in Fig. 1(a). It is given by the J ¼ 5=2; m J ¼ ± 3=2 j idoublet, where the quantization axis z is chosen as the local [111] axis of the cubic lattice 9 . The transverse components of the angular momentum thus have vanishing expectation values in the ground state, which in turn implies that they are invisible to the spin-flip scattering in the neutron-scattering experiments because Δm J = 3 excitations do not couple, to leading order, to the dipolar magnetic moment of neutrons. Using an effective pseudospin 1/2 representation of the ground state doublet 15 , one of the three components (historically denoted s y ) represents the octupolar moment, and the other two components (s x and s z ) transform as the the familiar dipole.
The key to the unusual properties of Ce 2 Zr 2 O 7 is the effective interactions between spin components, both dipolar and octupolar, belonging to the nearest-neighbors Ce ions, which are different from the thoroughly studied dipolar spin ice. A model nearest neighbor (nn) Hamiltonian incorporating all symmetryallowed spin-spin interactions on the pyrochlore lattice is 15 : with s i again expressed in the local frame, i.e., relative to the local [111] direction on a given Ce site (for a detailed expression in the global frame see Supplementary Note 1). J x , J y , J z , and J xz are magnetic coupling strengths and 〈ij〉 refers to nearest neighbor pairs of sites. While the form of the Hamiltonian is dictated by symmetry, the values of the parameters, which dictate the location of the material on a theoretical phase diagram, depend on its chemistry. Owing to the multiple competing scales, predicting these numbers accurately from first-principles calculations is known to be challenging and has been the subject of various past and ongoing debates (see for example, refs. 18,19 ). The route we take is to optimize the parameters in Eq. (1) by quantitatively reproducing results of multiple experiments using a combination of quantum and classical methods. After accomplishing this objective, we appeal to known theoretical predictions and infer the kind of phase Ce 2 Zr 2 O 7 realizes. We find that all our optimized parameter sets lie in a π-flux U(1) QSL part of a previously constructed mean-field phase diagram 20 .

RESULTS
Effective Hamiltonian from specific heat and magnetization By fitting the experimental magnetization and specific heat to (quantum) finite temperature Lanczos method (FTLM) calculations 21 (see Fig. 2, Methods and Supplementary Note 2 for more details), we have determined the parameters of the model Hamiltonian in Eq. (1). A crucial result from the modeling perspective is that we have identified the J y interactionwhich acts between the octupolar componentsto be the largest term (J y ≈ 0.1 meV), with two interactions (J x and J z ) playing a subleading role, and a vanishing J xz . We have obtained several sets of parameter values (see Supplementary Table 1) within the fitting error bars, shown by the black dots in Fig. 1(b), clustered around J y = 0.08 ± 0.01 meV, J x = 0.05 ± 0.02 meV, J z = 0.02 ± 0.01 meV. We show additional cost function analyses for a wide range of J x , J y and J z in Supplementary Note 3, that illustrate constraints on our fits given the current availability of experimental data. Additionally, since the fitting was performed using a 16-site cluster to reduce the computational cost, we also verified that the parameters adequately describe a 32-site system (see Supplementary Note 4).

Inelastic neutron spin structure factor
Note that since the octupolar s y moments do not couple to neutron spins to leading order, it is crucial to perform fits to the magnetization and specific heat as described above; attempts to fit solely the dynamical structure factors (see expressions in Supplementary Note 5) measured in inelastic neutron scattering (INS) are less reliable. We do employ the INS data, however, at the second stage of our fitting process. After determining H nn , we introduce an additional weak next-nearest neighbor (nnn) coupling J nnn (see Methods and Supplementary Note 1), which likely originates from the magnetic dipole-dipole interaction between Ce ions. Its optimal value J nnn~− 0.005J y ≈ −0.5 μeV was determined by comparing the Self-Consistent Gaussian Approximation (SCGA) prediction of the spin structure factor with INS data (see Fig. 2(c)).
Armed with this complete Hamiltonian (nn and nnn terms), we compute the energy-integrated and energy-resolved momentum-dependent neutron spin structure factor using classical Molecular Dynamics (MD) [22][23][24] for N = 1024 and N = 8192 sites. Representative comparisons with previous experiments on Ce 2 Zr 2 O 7 9,10 are shown in Fig. 1(c-e). Figure 1(c) shows the characteristic ring-like structure centered around Γ 0 = (0, 0, 0) point, with pronounced maxima at |q| = 1 (in particular q = (0, 0, 1) and ð 1 ffiffi 2 p ; 1 ffiffi 2 p ; 0Þ where we use the standard Miller indices (h, k, l) to denote the location in reciprocal space), consistent with experimental observations. Using the quantum-classical correspondence to rescale the MD data 23 , in Fig. 1(d) we show the results for a one-dimensional cross section in momentum space (Γ 0 → X ≡ (0, 0, 1) → Γ 1 ≡ (0, 0, 2)). The increased intensity at Fig. 1 Quantum spin liquid in a pyrochlore magnet and its experimental ramifications. a Depiction of a pyrochlore lattice with octupolar components forming 2-in-2-out ice states, and showing representative nearest neighbor (nn) and next-nearest neighbor (nnn) bonds. b Our model parameter sets are shown with black dots, superposed on the mean-field phase diagram of ref. 20 . c Dynamical structure factors at T = 0.06 K integrated over the energy range 0.00 to 0.15 meV obtained from classical molecular dynamics (MD) for 8192 spins along with the corresponding comparison with ref. 9 (T = 0.035 K) and ref. 10 (T = 0.06 K, with an additional background subtraction). The plots reported in the published references were adapted for the purpose of comparison with our results. Special points in the Brillouin zone are also shown. d, e The dynamical structure factor as function of energy and momentum using the MD data (for 1024 sites) compared to experiment. The quantum-classical correspondence dictates that the MD data be rescaled 23 by the factor βE (β = 1/k B T where k B is the Boltzmann factor, and E is the neutron energy transfer). Panel d is for the momentum (0, 0, l) cross section (Γ 0 → X → Γ 1 ) and panel e is the powder average. For panels c, d, e the numerically optimized Hamiltonian (parameter set no. 2 with both nearest neighbor and next-nearest neighbor terms, see Supplementary Table 1) was simulated. Additional Lorentzian convolution of width Γ was applied to mimic the limited experimental resolution, based on values reported in ref. 9 and ref. 10 . The color scheme used for the theoretical calculation and the two experiments differs, and in the absence of additional information only the relative variations should be compared. Note that the high intensity around q = (1, 1, ± 1) and (2, 2, − 2) seen in the experimental inset in panel c results from an imperfect subtraction (using high-temperature data) of the nuclear Bragg peaks, which are absent from our theoretical calculations. The same is true in panel e at |q|~1 Å −1 and |q|~2 Å −1 , where the maxima of intensity originate from imperfect subtraction of nuclear Bragg peaks.
the X point at low energies is broadly consistent with experimental findings (additionally, see Supplementary Note 6). Figure 1(e) shows our results for the powder averaged case, compared with the experimental data in ref. 10 , suggesting an overall agreement of the energy scales over the entire Brillouin zone (the X point corresponds to |q|~0.6 Å −1 in physical units, where the intensity is highest).
The fact that multiple experimental features are accurately reproduced by a model of dipolar-octupolar interactions between the Ce spin components in Eq. (1) (plus a small J nnn term) poses the question about what phase corresponds to its ground state. Indeed, the fact that the coupling J y between the nearest neighbor octupolar moments is antiferromagnetic and by far the largest suggests that the leading behavior is for the corresponding components to form a (classical) 2-in/2-out spin-ice manifold (In Supplementary Note 7, we show additional evidence for the importance of J y in explaining the experimental data). The presence of non-zero J x and J z interactions then adds quantum effects; generically, this opens the possibility of obtaining a quantum spinice phase.
Quantum spin-ice phase The phase diagram for the model Hamiltonian in Eq. (1) was studied in refs. 20,25 , using a combination of analytical and meanfield analysis as well as exact diagonalization. In our modeling of the magnetization and specific heat, we have determined four candidate sets of fitting parameters depicted by black dots in Fig. 1(b). All four fall deep into the parameter regime of the π-flux octupolar quantum spin-ice (π-OQSI) phase, based on the phase diagram of ref. 20 . In this phase, the emergent gauge field A takes a non-trivial ground state configuration that hosts a flux Φ = π through each hexagon [26][27][28] , where r i is the site of the dual diamond lattice, and the spins and A live on the links r i , r i+1 of the pyrochlore lattice, as shown schematically in the inset of Fig. 3(a). This fact follows from the effective U(1) quantum field theory, described in the Supplementary Note 8, with the flux-dependent contribution in the form where J ring $ ðJ x þ J z Þ 3 =ð64J 2 y Þ is positive and thus favors flux Φ = π in each hexagon, resulting in the π-OQSI phase.
Like "conventional" quantum spin ice, π-OQSI has gapless photons and two types of gapped excitations (magnetic and electric charges), in close analogy to Maxwell electrodynamics. Their approximate energy scales are illustrated in Fig. 3(a). The first type of gapped excitations are spinons created in pairs by flipping the octupolar moment on a single site, costing energy OðJ y Þ. These are analogs of the magnetic monopoles in electrodynamics, and observable in the specific heat as a characteristic Schottky peak at energy~J y ≈ 1 K, as our FTLM calculations corroborate in Fig. 2(a). The second type of gapped excitations corresponds to the so-called visons (analogs of electric charges), which are sources of the gauge flux violating the condition in Eq. (2). Note the energy scale for exciting visons at J ring~0 .01 K is very low. One therefore expects to find thermally excited visons even at the base temperature of the experiment, so that their gap, if not closed by the quantum dynamics, will not be separately resolved. Rather, visons will strongly interact and mix with the emergent gapless photons, named in analogy to the photons familiar from Maxwell electrodynamics, due to the overlap of their energy scales [cf. Fig. 3(a)].
The energy and temperature scales for observing the photons would correspondingly be very low. More importantly, they would not directly couple to neutrons because of the octupolar nature of the π-OQSI. To illustrate this point, we have computed the pseudospin correlation functionsS αβ ðqÞ ¼ 1 N P i;j e iqÁðriÀrj Þ hs α ðr i Þs β ðr j Þi in Fig. 3(b) by SCGA (see Methods for details). The largest components of this matrix are the diagonal onesS xx ,S yy andS zz . It is the octupolarS yy ðqÞ component in the middle panel of Fig. 3(b) that displays the pinchpoints characteristic of the spin-ice 4,29 , and the photons' gapless dispersion will emanate from the q location of those pinch-points. Crucially however, the s y components of spin, at leading order, do not couple to the magnetic field or to the neutron moment (see Methods, especially Eq. (6)), meaning that the aforementioned pinch-points will not feature in the experiment. Instead, the (energy-integrated) neutron-scattering structure factor is expressed in terms of true magnetic moments m μ = ∑ λ g μλ s λ that contain the g-factors, whose g μy components are all zero (see Supplementary Note 1). As a consequence,S yy drops out of the neutron structure factor, computed in Figs. 1(c) and 2(c) using Monte Carlo and SCGA respectively. The main contributors to the neutron-scattering structure factor are theS xx andS zz channels convoluted with the neutron-coupling form factor in Eq. (4). As a result, instead of the pinch-points, a sixfold, three-rod-crossing-like pattern is observed. Such a rod pattern is expected for a pyrochlore lattice with nearest-neighbor interactions only 30 and is associated with the dispersion of spinons (magnetic monopoles) in the context of quantum spin ice [30][31][32] . Upon inclusion of (weak) nnn interactions J nnn , the pattern of crossing rods deforms into a characteristic ring-like structure that appears in S(q), as shown in the three panels of Fig. 2(c). Phenomenologically fixing the value J nnn~− 0.5 μeV matches very well with the experimental INS observations [see Fig. 1(c)]. Thence, while it is tempting to associate the intensity variation along the rods in Fig. 2(c) as the disappearance of pinch-point intensity centered at (1, 1, 1) (and equivalent points), which has been predicted to be the quintessential feature of the dispersive, quantum photon modes in dipolar quantum spin ice 33 , our analysis rather suggests that these are the consequence of the small nnn interactions between Ce ions that modulate the INS intensity along the [111] direction. As for the emergent photons, while they are indeed expected to be present in the π-OQSI phase, as shown in Fig. 3(a), their octupolar nature turns out to render them much less visible to neutrons, and they can only be detected via weaker, higher-order coupling to neutrons at large momentum transfer 34,35 , or indirectly, for instance through their contribution to the low-temperature specific heat (at T ≲ J ring /k B ≈ 0.01 K).

DISCUSSION
Given the effective microscopic description of Ce 2 Zr 2 O 7 , an interesting question is to ask how the (putative) octupolar quantum spin-ice state responds to the application of an external magnetic field 25 . While the octupolar s y pseudospin components do not couple linearly to the field and will remain in the 2-in/2-out configuration, the s z components will cant along the field direction (see Eq. (6), where g x ≈ 0). In order to elucidate the experimental consequences, we compute the in-field spin structure factor S(q) within classical Monte Carlo calculations on our model for 8192 sites, shown in Fig. 4. As a function of increasing field along the [001] direction, the ring-like structure in S(q) quickly weakens (Fig. 4b), until eventually disappearing and giving way to sharp Bragg peaks in high fields (Fig. 4c). These predictions are to be compared with future INS data in an applied magnetic field.
The fact that the magnetic octupolar degrees of freedom do not couple in the leading order to neutron spin or to the external magnetic field, makes the octupolar spin liquid difficult to detect. Its elusive nature may however prove to be a blessing in disguise, as a reduced coupling to magnetic defects and associated stray magnetic fields-which are known to destabilize the more conventional dipolar spin liquids-is similarly suppressed. Indeed, chemical disorder on magnetic sites is believed to be the leading reason for the failure to observe the quantum spin liquid behavior in, for instance, the herberthsmithite kagome compounds despite their high crystallographic quality 36 . The fact that neither magnetic order nor spin glassiness is seen in cerium pyrochlores Ce 2 Zr 2 O 7 and Ce 2 Sn 2 O 7 34,37 may be taken as an additional, albeit indirect, evidence of the robustness of the underlying octupolar spin liquid. While the possibility of such a quantum spin liquid has been entertained in seminal theoretical studies before 15,20 , our present work firmly identifies Ce 2 Zr 2 O 7 as a very promising host for the π-flux octupolar quantum spin-ice phase. The present study also underscores the importance of carefully fitting multiple experiments, including specific heat and magnetization, in addition to the INS spectra, to determine the effective model Hamiltonian, which otherwise may be plagued with uncertainties affecting the searches and identification of QSLs 18,19 .

Fitting model Hamiltonian parameters
In addition to the nearest neighbor (nn) Hamiltonian in Eq. (1) we have considered the next-nearest neighbor (nnn) interaction, which in the local basis is given by, where 〈〈i, j〉〉 refers to nnn pairs of sites. For the case of an applied external magnetic field, the Zeeman term must also be accounted for. This term involves the coupling of magnetic field to effective spin 1/2 degrees of freedom, which are not the usual and familiar dipoles, and are instead dipolar-octupolar doublets. A key observation is that the octupolar magnetic moments do not couple, to linear order, to the external magnetic field. The octupolar moments can in principle couple to the third power of the magnetic field H 3 (ref. 20 ), however, this effect is expected to be small for the experimentally accessible fields. It can be shown that only the component of the field H along the local z-axis couples to the dipole moment 15 , as follows: where |h| = μ B μ 0 H is the effective magnetic field strength, μ B is the Bohr magneton. Note that when projected onto the J = 5/2 multiplet, one expects the Landé g-factor g z = 2.57 and g x = 0. However, if an admixture of higher spin-orbit multiplet (J = 7/2) is present in the ground state doublet, one generically expects a non-zero value of g x (see Supplementary Note 9), which we allow for in our modeling. As explained in the main text, the Hamiltonian parameters in Eq. (1) were determined using both zero and finite applied magnetic field data. For the specific heat and magnetization, we have performed quantum FTLM calculations on a 16-site cluster. Details of the general technique can be found in ref. 21 and previous application to some pyrochlore systems can be found in ref. 38 . Our convergence checks of the method (see Supplementary Note 2) suggest that the statistical errors present in the specific heat and magnetization are smaller than the line width used in Fig. 2. Figure 2(a) shows the temperature dependence of the specific heat and suggests that fitting it is somewhat challenging, both in zero and applied field. Our four distinct parameter sets are generally able to describe specific heat very well at higher temperature T ≳ 1 K, however, the lower-temperature behavior is trickier and none of the parameter sets used satisfactorily fit the data especially in the absence of applied magnetic field. We attribute this to the finite-size effects in our FTLM calculations, which become more pronounced at lower temperatures. The magnetization as a function of applied field strength, shown in Fig. 2(b), appears to be less prone to finitesize effects (also see Supplementary Note 4) and is fitted reasonably well, especially at low fields. The discrepancy at high fields is not entirely unexpected, previous experimental reports have suggested a changing g-factor past a μ 0 H = 4 T field 9 , an effect not built into our model.
Despite the finite-system size limitations of the numerical calculations, and the finite energy and momentum resolution in INS experiment, all our parameter sets obtained are in agreement with J y being antiferromagnetic. For three of these sets, we find J y is large compared to the other two interactions (J x and J z ) in Eq. (1), as depicted in Fig. 1(b) (also see Supplementary Note 7). A fourth parameter set, obtained by constraining g x to be zero, suggests J y ≈ J x > J z . To further build confidence in our results, we have performed a brute force scan of J x and J z for representative fixed values of J y and constructed a contour map of an appropriately defined cost function (see Supplementary Note 3, more subtleties with the fitting and additional competitive parameter sets are also discussed). Additional future experiments could potentially further constrain the values of these coupling constants.
The second step of our parameter fitting involved the determination of J nnn , which we found to be small relative to J y . Despite its smallness, it is responsible for significant reorganization of intensity in the Brillouin zone. Figure 2(c) shows the static structure factor in the (h, h, l) plane computed with SCGA (the details of which will be explained shortly) for representative values of J nnn . Performing a brute force line search, we determined the optimal value J nnn = − 5 × 10 −4 meV working in steps of 10 −4 meV.
The following values of the parameters (which we refer to as set no.

Details of the Monte Carlo and molecular (spin) dynamics calculations
For the Hamiltonian H H nn þ H nnn þ H Z , the dynamical structure factor is computed by integrating the classical Landau-Lifshitz equations of motion which describes the precession of the spin in the local exchange field. We carry out all our calculations by transforming our Hamiltonian to the global basis, and have used the label S i to represent spins in this basis (see Supplementary Note 1 for more details on transformations between local and global bases).
Following the protocol adopted in previous work 22,23 , the initial configuration (IC) of spins is drawn by a Monte Carlo (MC) run from the Boltzmann distribution expðÀβHÞ at temperature T = β −1 = 0.06 K. Then, for each starting configuration the spins are deterministically evolved according to Eq. (7) with the fourth order Runge-Kutta method. This procedure, referred to as molecular dynamics (MD) (also see Supplementary Note 5), is repeated for many independent IC (their total number being N IC ) and the result is averaged, We perform a Fourier transform in spatial and time coordinates to get the desired dynamical structure factor. We work with N = 16L 3 sites, where L 3 is the number of cubic unit cells, the results in Fig. 1(c) are for L = 8 (N = 8192) and L = 4 (N = 1024) for Fig. 1(d, e). N IC ≈ 10 4 was used, and each IC was evolved for 500 meV −1 in steps of δt = 0.02 meV −1 . The statistical error bars in the classical MD/MC simulations presented are small enough to not change the features of the color plots.
To obtain an estimate of the quantum dynamical structure factor, we used a classical-quantum correspondence, which translates to a simple rescaling of the classical MD data by βE. More details and justification can be found in ref. 23 .
Convolution with Lorentzian function to mimic limitations of experimental resolution Figure 1(d, e) shows significant broadening along the energy axis, an effect not captured to the same extent by the raw (rescaled) MD data. Since the interaction energy scales in the material are small, a fairer comparison between experiment and theory is achieved by modeling the instrument's energy resolution. Using Γ values in the ballpark suggested by refs. 9,10 , we convolved our data with a Lorentzian factor, The integral was approximated by a sum over discrete E 0 points in steps of 0.01 meV.

Details of the SCGA calculations
The Self-Consistent Gaussian Approximation is an analytical method that treats the spin in the Large-N (flavor) limit. Our calculation follows closely ref. 4 . In this study, we first treat s x,y,z as independent, freely fluctuating degrees of freedom. The Hamiltonian in momentum space is written as where S ¼ ðs x 1 ; s x 2 ; s x 3 ; s x 4 ; ; s z 3 ; s z 4 Þ: The interaction matrix H Large-N is the Fourier transformed interaction matrix that includes the nearest and nextnearest neighbor interactions.
We then introduce a Lagrangian multiplier with coefficient μ to the partition function to get in order to impose an additional constraint of averaged spin-norm being one, or hs 2 1 þ s 2 2 þ s 2 3 þ s 2 4 i ¼ 1: For a given temperature k B T = 1/β, the value of μ is fixed by this constraint via relation Z BZ d k where λ i (k), i = 1, 2, … , 12 are the twelve eigenvalues of βH LargeÀN . With μ fixed, the partition function is completely determined for a free theory of S, and all correlation functions can be computed from βH LargeÀN þ μI Â Ã À1 .

DATA AVAILABILITY
The data analyzed in the present study is available from the first author (A.B.) upon reasonable request.