A pseudofermion functional renormalization group study of dipolar-octupolar pyrochlore magnets

Motivated by recent experiments on Ce$_2$Zr$_2$O$_7$ that reveal a dynamic, liquid-like ground state, we study the nearest neighbor XYZ Hamiltonian of dipolar-octupolar pyrochlore magnets with the pseudofermion functional renormalization group (PFFRG), which is numerically implemented by the SpinParser software. Taking the interaction between the octupolar components to be dominant and antiferromagnetic, we map out the phase diagram demarcating the quantum disordered and magnetically ordered states. We identify four distinct phases, namely the $0$-flux and $\pi$-flux quantum spin ices, and the all-in-all-out magnetic orders along the local $z$ and $x$ axes. We further use the static two-spin correlations output by the PFFRG algorithm to compute the polarized neutron scattering cross-sections, which are able to capture several qualitative features observed experimentally, in the materially relevant parameter regime that stabilizes the $\pi$-flux quantum spin ice. Our results provide support for a quantum spin liquid ground state in Ce$_2$Zr$_2$O$_7$.


I. INTRODUCTION
A spin liquid is, roughly speaking, a cooperative paramagnet with disordered yet correlated spins.Depending on whether its dynamics are governed by thermal or quantum fluctuations, a spin liquid is primarily classified as classical or quantum.A prominent example of spin liquids is spin ice [1] on the threedimensional pyrochlore lattice, where each unit tetrahedron displays a two-in-two-out spin configuration.While there exist concrete experimental evidences of classical spin ice [2][3][4] in Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 [5][6][7][8][9], the detection and confirmation of quantum spin ice (QSI) [10][11][12][13] in candidate materials remains an ongoing work.QSI is of fundamental theoretical interest as it realizes a lattice analog of quantum electrodynamics, which exhibits gapless photons and two types of gapped excitations, namely spinons and visons, which are the electric and magnetic monopoles of the emergent  (1) gauge theory.
In this paper, we theoretically investigate the  = 1/2 XYZ model of the dipolar-octupolar pyrochlore magnets using the pseudofermion functional renormalization group (PFFRG) , which is numerically implemented by the SpinParser software [60,61].As the name itself suggests, the two essential components of PFFRG are (i) a pseudofermion representation of the spin operator, so that a Hamiltonian with arbitrary twospin interactions is cast into an interacting fermion problem, followed by (ii) a functional renormalization group [62][63][64] analysis.The purpose is to obtain a low-energy description of the system in terms of the renormalized self-energy and twoparticle vertex function, which are used to compute the static component of the magnetic susceptibility.Compared to other methods, PFFRG has the advantage of treating the thermodynamic limit by effectively truncating the interaction range of fermions, and it is free of the sign problem encountered in quantum Monte Carlo.Therefore, PFFRG is advocated as a suitable tool to study frustrated magnetism in three spatial dimensions.
We first map out a phase diagram in the parameter space relevant to Ce 2 Zr 2 O 7 , as shown in Fig. 1.Tracking the evolution of the static susceptibility with respect to the RG cutoff, PFFRG is able to distinguish symmetry-breaking ordered states and symmetry-preserving paramagnetic states [51,59,60], the latter of which are putative quantum spin liquids at low temperatures.To further differentiate among the quantum spin liquids, we look for qualitative differences in the momentum-resolved static susceptibilities across the parameter space, and use insights gained from previous theoretical analyses [30,33,35].Within the parameter space of interest, our PFFRG analysis identifies four distinct phases: a 0-flux QSI, a -flux QSI, and two all-in-all-out (AIAO) magnetic orders along the local  and  axes, respectively [65].The proposed parametrizations of Ce 2 Zr 2 O 7 in the existing literature [16,18,36] are found to lie within the -flux QSI phase.
Although the static susceptibility is not exactly the equaltime spin structure factor, the former is a good estimate of the latter if the spectral weight of the dynamical spin struc-  [16,36] at the level of nearest neighbor interactions.The black dots demarcate the symmetry-preserving and symmetrybreaking regions as signalled by the breakdown of smooth renormalization group flow, with the error bars mainly reflecting the resolution of the grid for the calculations.The lines (without dots) separating the two spin ices and the two magnetic orders are obtained by further examining the diagonal components of the static susceptibility.See Sec.IV A for more details.ture factor is concentrated at low energies, or the dynamical spin structure factor is nonzero only within a relatively narrow range of finite energies.In principle, one can calculate the magnetic susceptibility at arbitrary (Matsubara) frequency, but the integration over frequency to obtain the equal-time spin structure factor leads to additional numerical errors [57].As an approximation, we thus use the static susceptibility itself to compute the neutron scattering cross-section.We back this approximation with calculations of the static susceptibility and the equal-time spin structure factor using gauge mean field theory [32,33], which reveal highly similar intensity profiles.We find that the PFFRG calculated neutron scatterings at several parameters stabilizing the -flux QSI not only reproduce the rod motifs in the spin-flip channel, but also capture the positions of the intensity minima in the non-spin-flip channel, as seen in the experimental data reported by Ref. [16].Our results thus offer support to the case of a quantum spin liquid ground state in Ce 2 Zr 2 O 7 , which may be the -flux QSI.
The rest of this paper is organized as follows.Sec.II describes the XYZ model of dipolar-octupolar pyrochlore magnets.Sec.III briefly explains the PFFRG methodology, in particular how to distinguish ordered and disordered states by analyzing the cutoff dependence of the static susceptibility.Sec.IV presents the results from the PFFRG analysis, namely the different phases and the calculated structure factors.Sec.V summarizes our work, discusses our results in light of other existing works, and comments on possible future directions.

II. MODEL
Dipolar-octupolar pyrochlore magnets such as Ce 2 Zr 2 O 7 consist of interacting pseudospin-1/2 degrees of freedom, of which the  ( and ) components transform as magnetic octupoles (dipoles) under lattice symmetries and time reversal [26], on a three-dimensional network of corner-sharing tetrahedra.The most generic nearest neighbor Hamiltonian is given by [26] where the pseudospin-1/2 components are defined according to the local coordinates at each site, in which the  axis points from the center of a tetrahedron to one of its vertices.The offdiagonal term   in (1) can be removed by a rotation about the  axis, which results in the XYZ model, We will drop the tildes in (2) as it is suggested that   ≈ 0 in the candidate material [16,18,36] The pyrochlore lattice is a face-centered cubic (fcc) lattice with four sites per unit cell.We use the following convention for the primitive translation vectors and the sublattice displacements where x, ŷ, ẑ are three orthonormal vectors in the global frame.

III. METHOD
We employ the pseudofermion functional renormalization group (PFFRG)  to study the  = 1/2 XYZ Hamiltonian (2) on the pyrochlore lattice.We only provide a brief description of the method in this section, with some further details given in Appendix A. Interested readers may refer to, e.g., Refs.[39,51,[57][58][59][60], for an in-depth discussion.

A. Pseudofermion Functional Renormalization Group
We start by representing the spins in terms of pseudofermions, where Eq. ( 5) is a faithful representation only when the single occupancy constraint  † ↑  ↑ +  † ↓  ↓ = 1 is satisfied locally.To a good approximation, this constraint is enforced on average by setting the chemical potential to be zero, as its violation leads to energetically unfavorable  = 0 local defects, which are thermally suppressed at low temperatures [53,59,60].
The strongly interacting pseudofermion Hamiltonian ( 6) is then subject to the functional renormalization group (FRG) analysis [62][63][64].The central objects of FRG are one-line irreducible vertex functions, or simply vertices, which encode the effective -particle interactions [57,66].The FRG flow equations are generated by introducing an infrared cutoff Λ in the Matsubara frequency to the bare propagator  0 () = 1/, such that With the magnetic couplings of the original spin model treated as bare interactions at Λ −→ ∞, we are ultimately interested in the low-energy effective theory at Λ −→ 0. Differentiating the generating functional of the vertices with respect to the cutoff, one obtains an infinite hierarchy of coupled integro-differential equations, in which the flow of the -particle vertex involves vertices up to the ( + 1)-th order.We then apply the Katanin truncation scheme [67], so that we only have to solve the flow equations for the one-particle vertex, which is equal to the self-energy up to a minus sign, and the two-particle vertex while partially retaining the effect of the three-particle vertex [58,59].The structures of these vertices are considerably simplified by the symmetry of the original spin model as well as the gauge redundancy from the pseudofermion construction [51,59].For instance, the self-energy is an imaginary and antisymmetric function that depends only on the Matsubara frequency.More details can be found in Appendix A.

B. Numerical Implementation
All PFFRG calculations in this work are performed with the newly available SpinParser software [60,61], which solves the flow equations numerically.As shown in Appendix A, the flow equations involve summations over Matsubara frequencies and lattice sites, see (A5a) and (A5b).The Matsubara frequency becomes continuous in the zero temperature limit, while the number of lattice sites grows to infinity in the thermodynamic limit.Further approximations are thus required for a numerical solution.
To this end, the frequency axis is rediscretized such that frequency dependent quantities are evaluated for a finite set of frequencies, supplemented with a linear interpolation scheme.We choose   = 144 frequencies distributed logarithmically around  = 0, with || max = 100 and || min = 0.001.The numerical integration over frequency is performed with a trapezoidal scheme [60].On the other hand, we set the two-particle vertex to be zero if the two lattice sites involved are further apart than  = 6 nearest neighbor bonds.Such a finite truncation range allows us to effectively study the thermodynamic limit without imposing specific boundary conditions.Finally, the cutoff is decreased in steps via Λ +1 = Λ  with the factor  < 1.We choose  = 0.98, and the initial (final) cutoff at Λ  = 100 (Λ  = 0.01), which is much greater (smaller) than any intrinsic energy scale.All these specifications can be done within SpinParser [60,61].

C. Magnetic Susceptibility
A useful physical observable that can be calculated from PFFRG is the static component ( = 0) of the magnetic susceptibility (two-spin correlator), We emphasize that the magnetic susceptibility is also a function of the cutoff Λ, though it is not shown explicitly.
Analyzing the Fourier-transformed static susceptibility allows us to infer the ground state of the system as follows.The PFFRG calculation assumes the full symmetry of the Hamiltonian, while a magnetically ordered state corresponds to spontaneous symmetry breaking.The onset of magnetic order causes a breakdown of the FRG flow, which typically manifests as a divergence or a kink in the cutoff dependence of the magnetic susceptibility [57,60].To determine the ground state, one can trace the evolution of (k) ≡    (k) as Λ decreases at some momentum k = k * , which is typically chosen such that (k * ) is largest.If (k * ) becomes nonanalytic at some critical cutoff Λ c , then it signals a transition into a symmetry broken phase, and k * right before the flow breakdown is taken as the ordering wavevector.In contrast, if (k * ) remains smooth and finite down to Λ −→ 0, then it indicates a paramagnetic ground state that preserves all symmetries.It is worth remarking that the solutions of the flow equations for Λ < Λ c are no longer physically meaningful due to symmetry breaking, so one would not be able to access the true ground state at Λ −→ 0 should there be a multistep ordering process [53].Some degree of uncertainty is inevitable in locating the critical cutoff Λ c by inspection.Moreover, it may be difficult to detect the nonanalyticity in (k * ) when the phase transition takes place at some small cutoff value Λ ≈ 0. To complement the procedure described in the previous paragraph, we compare the on-site susceptibility   ≡     ( = 0) as a function of Λ for multiple truncation ranges  [52,53].The rationale is that paramagnets and spin liquids only exhibit short-range correlations, so that   converges already at small .In contrast, long-range correlations become important when there is a tendency for magnetic ordering, which results in an unambiguous discrepancy between   for small and large .One can then set a quantitative threshold for the discrepancy to define Λ c [53].

A. Phase Diagram
We have outlined in Sec.III C how magnetic orders are distinguished from paramagnetic phases such as quantum spin liquids in general PFFRG calculations.We now specialize to the pyrochlore XYZ model (2), and discuss how to further differentiate among magnetic orders, as well as quantum spin liquids.
When a symmetry breaking phase transition occurs, k = k * that yields the maximum of (k) at Λ = Λ c is taken as the ordering wavevector.Plotting (k) at the critical cutoff over an extended region in the reciprocal space, one typically observes a rather sharp peak at k * , see Fig. 3l for instance.Different ordering wavevectors correspond to distinct magnetic orders and thus serve as primary labels of the latter.
However, in the parameter space of interest, we find only k = 0 orders, which signify ferromagnetic correlations, when   or   is sufficiently negative.A further distinction is made by examining whether the  correlation at the Γ point dominates over , or vice versa.Along the line   =   , we have    (k) =   (k) for all k, see Fig. 2b  of a tetrahedron are either all aligned or all antialigned to their respective local  axes, which is known as the all-in-allout (AIAO) order.This state is labeled as Z-AIAO following the convention in Refs.[35,37], to distinguish it from the X-AIAO state where    (k = 0) is dominant.The Z-AIAO and X-AIAO magnetic orders are naturally separated by the   =   line, see Fig. 1.
It is much more challenging to distinguish quantum spin liquids, as they preserve the full symmetry of the Hamiltonian.While PFFRG is able to identify spin liquid ground states, it does not by itself offer a classification of different spin liquids.Nevertheless, we can try to analyze the static susceptibilities calculated from PFFRG and see if there exists any qualitative difference among them.We also make use of insights developed in other theoretical studies, which identify the 0-flux and -flux quantum spin ices (QSIs) as the two major spin liquid candidates in the parameter region of interest [30,33,35].In the perturbative regime where   ≫ |  |, |  |, it has been well established [10] that the 0-flux (-flux) QSI are stabilized for Hence, we plot and examine the diagonal components of   (k) in the [ℎℎ] plane in reciprocal space.Momenta are measured in the reciprocal lattice units; e.g., (ℎ, ℎ, ) = (1, 1, 3) means k = 2( x + ŷ + 3ẑ)/.For every parameter that stabilizes a quantum spin liquid ground state,    (k) shows bowtie-like motifs signifying spin ice correlations, with pinch point singularities at (0,0,2) and (1,1,1) that are broadened, see Figs. 3b and 3f.Meanwhile,    (k) and   (k) in the spin liquid regime largely fall into two categories: both of them either show (i) diffuse peaks at the Γ point, see Figs. 3a and 3c, or (ii) the bowtie patterns similar to those in   (k) but without obvious pinch points, see Figs. 3e and 3g.We associate (i) and (ii) with the 0-flux and -flux QSIs, respectively [33,37].
These two QSIs are separated approximately by the line   +   = 0 (see Fig. 1), which is consistent with perturbation theory.While there is a rather thin slice of the 0-flux QSI in the phase diagram, the -flux QSI occupies a much larger area, likely owing to stronger frustration from positive transverse couplings.Within the -flux QSI, the bowtie motifs in the  correlation are sharper (more diffuse) for larger (smaller)   , and the same is true for the  correlation and   , when the color scale is chosen to range from zero to the maximum intensity.For example, one can notice the difference between Figs. 3g (  = 0.1) and 5d (  = 0.7).
We remark that there exists a small parameter region in which    (k) displays a maximum at the Γ point while   (k) displays the bowtie patterns, or vice versa, see Figs. 8a, 8c, 9e, and 9g in Appendix B. As this only appears near the phase boundary between the -flux QSI and the 0-flux QSI or one of the magnetic orders, we interpret it as a tendency of the -flux QSI to develop ferromagnetic correlations in proximity to a phase transition, and refrain from attributing it to a new kind of quantum spin liquid.Finally, we caution that it is possible for distinct quantum spin liquids to exhibit highly similar two-spin correlation functions that would not be resolved by PFFRG.
Throughout this paper, intensity plots in the [ℎℎ] plane are calculated at the smallest cutoff Λ = 0.01 for quantum spin liquids, and at the critical cutoffs Λ c for magnetically ordered states.

B. Neutron Scattering Cross Sections
To compute the energy-integrated neutron scattering cross section for comparisons with experiments, one should, strictly speaking, use the equal-time spin structure factor, which is formally given by the integral of the magnetic susceptibility (8) over the Mastubara frequency, However, this additional frequency integration, which is performed over a discrete mesh, leads to further numerical errors [57].Therefore, the existing PFFRG literature mostly focuses on analyzing the static susceptibility instead of the equal-time spin correlator.We will resort to the static susceptibility [69] and back our approximation with calculations using gauge mean field theory.At very low temperatures, the momentum-resolved static susceptibility ( 9) is related to the dynamical spin structure factor S  (k, ) via the Kramers-Kronig relation and the fluctuation-dissipation theorem [57,70,71], On the other hand, integrating the dynamical spin structure factor over energy, one obtains the momentum-resolved equaltime spin structure factor, Note that   (k) and S  (k) defined above measure the twospin correlations in the local coordinates.
To make connections with experiments, one should consider the (total) neutron scattering structure factor where ẑ is the unit vector along the local  axis at site , assuming that external magnetic fields and neutrons only couple to the local  components of the pseudospins.In polarized neutron scattering ( 13) is further decomposed into spin-flip (SF) and non-spin-flip (NSF) channels [12,16,37], where û is the unit vector along the direction of the neutron polarization, which is perpendicular to the scattering plane, and v(k) = ( û × k)/| û × k|.
Due to the numerical uncertainties in calculating S    as mentioned in the beginning of this subsection, we replace it by     in ( 13), (14a), and (14b), and denote the resulting total, spin-flip, and non-spin-flip neutron scattering structure factors by  TOT (k),  SF (k), and  NSF (k), respectively.We then compare our results to the experimental data reported in Ref. [16].While the static susceptibility is not quite the equal-time spin correlator, it is a reasonable estimate of the latter if the spectral weight of S(k, ) is concentrated at low energies, or S(k, ) is nonzero only within a relatively narrow range of finite energies, by (11) and (12).We further support the replacement by a direct comparison between the static susceptibility and the equal-time spin structure factor calculated from the gauge mean field theory in Sec.IV C.
We first examine the two proposed parametrizations for Ce 2 Zr 2 O 7 in the existing literature.Ref. [16] finds (  ,   ,   ) = (0.063, 0.064, 0.011) meV up to a permutation of   and   , which is supported by a subsequent work [18].On the other hand, Ref. [36] finds a number of parametrizations clustered in the region (  ,   ,   ) = (0.05 ± 0.02, 0.08 ± 0.01, 0.02 ± 0.01) meV.Scaling   to 1, we get (  ,   ) ≈ (1, 0.17) and (0.63, 0.25), both of which lie within the -flux QSI phase, see Fig. 1.The calculated spin-flip and non-spin-flip neutron scatterings are shown in Figs.4b, 4c, 4h, and 4i.We find that the intensity variation of  NSF (k) is generally confined to a narrow range, so the profile of  SF (k) highly resembles that of  TOT (k).The latter is thus not shown separately.
The most eminent feature in  SF (k) is the rod-like distribution of high intensities [72], which is seen in the polarized neutron scattering experiment [16] as well as reproduced in a number of theoretical calculations [16,33,36,37].The rods also show apparent narrowing and widening in certain sections while remaining connected, see Figs. 4b and 4h.In particular, the vertical rod is narrowed at (0, 0, 2) and widened at (0, 0, 3), resembling a "neck" and a "head".The experimental data in Ref. [16] also features a "neck" at (0, 0, 2), albeit much narrower like a pinch point.Moreover, the narrowing and widening of the rods in other directions are less severe, which is consistent with the experiment.
We now turn to  NSF (k).The intersection between the Brillouin zone (BZ) boundaries of the fcc lattice and the [ℎℎ] plane is a network of edge-sharing hexagons, which looks like a honeycomb lattice compressed along one of the bond directions.We observe that the minima of  NSF (k) are located at the centers of these hexagons, see Figs. 4c and 4i, which is consistent with the experiment.However, the maxima do not take place exactly along the BZ boundaries as seen in the experiment, but rather form a stripe-like pattern.We find that increasing   leads to a starker contrast between the narrowed and widened sections of the rods in  SF (k), as well as a more diffuse intensity background for the minima in  NSF (k), see Figs. 4e, 4f, 4k, and 4l for calculations at the parameters (  ,   ) = (1, 0.5) and (0.63, 0.37).These changes coincide with the sharpening of the bowtie motifs in   (k): tighter "knot" of the bowtie is accompanied by narrower "neck" of the rods, c.f. Figs.4a and 4d.While introducing further nearest neighbor interactions may substantially improve the agreement between theory and experiment, as demonstrated by Ref. [36], we have shown that the nearest neighbor XYZ model is able to qualitatively capture the main features of the polarized neutron scattering cross-sections reported in Ref. [16].
Additional calculations of the spin-flip and non-spin-flip neutron scatterings at other parameters reveal that the -flux QSI with   ≳ 0.3 generically displays both (i) high intensity rods with the head-and-neck features in the spin-flip channel, and (ii) well defined minima at the BZ centers in the non-spinflip channel, see Figs. 5b, 5c, 5e, and 5f for instance.For smaller or negative values of   , we can still observe the rod motifs in the spin-flip channel, but the intensity profile of the non-spin-flip channel becomes more stripe-like, which looks like Fig. 4c.
Finally, we show plots of  SF (k) and  NSF (k) calculated at two choices of parameters that stabilize the 0-flux QSI, see Figs. 5h, 5i, 5k, and 5l.Roughly speaking, the intensity distribution of the 0-flux QSI is opposite to that of the -flux QSI: (i) the rods now carry low instead of high intensities in  SF (k), while (ii) the BZ centers are now maxima instead of minima in  NSF (k).These distinctions, which are also observed in two other theoretical studies [33,37], can serve as criteria in differentiating the 0-flux and -flux QSIs in future neutron scattering experiments.

C. Gauge Mean Field Theory
When computing the neutron scattering cross sections in Sec.IV B, we have replaced the equal-time spin correlator with the static susceptibility, as PFFRG is unable to calculate the former accurately.Here, we support the validity of this replacement by calculating and comparing the static susceptibility and the equal-time spin structure factor using gauge mean field theory (GMFT) [32,33].We can further compare the GMFT results in this subsection with the PFFRG results in the previous subsection.Both show the bowtie patterns in   (k), the rod-like distributions of high intensities with the head-and-neck features in  SF (k), and the minima at the BZ centers in  NSF (k).In addition, GMFT reveals (i) sharp point-like maxima along the highintensity rods in the spin-flip channel, so that an apparent dent of intensity is seen in the vicinity of k = 0, as well as (ii) maxima along the BZ boundaries in the non-spin-flip channel.The agreement between the GMFT results and the experimental data in Ref. [16] is excellent.

V. DISCUSSION
In summary, we have employed the pseudofermion functional renormalization group (PFFRG) to study the nearest neighbor XYZ model on the pyrochlore lattice, in the parameter regime relevant to the quantum spin liquid candidate Ce 2 Zr 2 O 7 .PFFRG analyses of pyrochlore magnets in the existing literature have mostly focused on Heisenberg models [50,54,[73][74][75], which are isotropic in spin space.Applications of PFFRG to anisotropic spin models on the pyrochlore lattice only appeared quite recently; our work contributes to this effort in addition to two others that investigate the Heisenberg-Dzyaloshinskii-Moriya model [76] and the non-Kramers pyrochlore model [77].
We present a phase diagram that contains two quantum spin ices (QSIs) and two magnetically ordered states, and the phase boundaries largely agree with existing theoretical works [30,33,35].Approximating the equal-time spin structure factor by the static susceptibility, we compute the spin-flip and non-spin-flip channels of the neutron scattering cross-sections at various parameters.We back such an approximation with calculations using gauge mean field theory (GMFT) [32,33].We find that the computed neutron scattering cross-sections are able to reproduce several qualitative features seen in the experimental data of Ref. [16] across a wide range of parameters within the -flux QSI phase.In other words, we have demonstrated that a reasonable agreement with the neutron scattering experiment can already be obtained at the level of nearest neighbor interactions, though our results may be further refined by including second nearest neighbor interactions as proposed in Ref. [36].More importantly, our results support the case of a quantum spin liquid ground state in Ce 2 Zr 2 O 7 , which is likely the -flux QSI.
Apart from PFFRG, theoretical methods such as numerical linked cluster [16], molecular dynamics [37], exact diagonalization [37], and gauge mean field theory [33] are able to reproduce the rod motifs seen in the spin-flip channel.However, numerical linked cluster and molecular dynamics do not capture the intensity variation in the non-spin-flip neutron channel, while exact diagonalization and gauge mean field theory do.GMFT currently yields the best agreement with the polarized neutron scattering experiment [16] among these methods.A rather intriguing feature of the GMFT calculations is that the mean field amplitudes associated with the non-  -conserving terms converge to zero, which effectively reduces the XYZ model to an XYX model [33].Such an emergent  (1) symmetry is corroborated by exact diagonalization results [37], which find S   (k) = S  (k) even though   ≠   .On the other hand, molecular dynamics [37] and the PFFRG analysis in this work predict that S   (k) ≠ S  (k) and    (k) ≠   (k) for   ≠   .It will be interesting to resolve this disagreement in future investigations.
We also point out several other possible directions for future studies.As mentioned in the main text, the pseudofermion representation of spins in PFFRG introduces unphysical states with zero or double occupancies, which are argued to be thermally suppressed in the zero temperature limit.However, they may become important at finite temperatures, which renders the PFFRG application inaccurate.To overcome this problem, the pseudo-Majorana functional renormalization group (PM-FRG) [74,78,79] has recently been developed, where the spin operator is represented in terms of SO(3) Majorana fermions.This construction generates no unphysical states, but rather redundant physical states due to a Z 2 gauge freedom.PMFRG allows the calculations of thermodynamic quantities such as free energy and heat capacity as a function of temperature.It will be interesting to apply PMFRG to the pyrochlore XYZ model.For the truncation of the infinite hierarchy of FRG integro-differential equations, one may also apply the more elaborate multiloop scheme [75,80] rather than the Katanin scheme.
From a broader perspective, it will be desirable to obtain the dynamical spin structure factor directly from PFFRG.This is not possible with the current PFFRG scheme, which is formulated in imaginary time, as the analytic continuation from Matsubara to real frequencies is a challenging numerical problem [58,59,81].The ability to calculate the dynamical spin structure factor would make PFFRG a more powerful theoretical tool in the study of frustrated magnetism, given the importance of inelastic neutron scattering experiments.ceptibility to study quantum spin liquids at  −→ 0 here, it is sensible to look nonetheless at the width of the pinch point as a proximate measure of the inverse correlation length and infer the typical separation between excitations [50].

𝜒
For concreteness, we study the XY and XXZ models with  ∈ [0, /2].For each  = , , , we calculate   (k, Λ = Λ  ) along a one-dimensional cut k = (0, 0, ) through k = 2X = (0, 0, 2), at which the pinch point is centered.To approximate the corresponding correlation length   , we first extract the half width at half maximum (HWHM  ) of the intensity along this cut [50], and plot it as a function of , see Figs. 11a and 11b.We note that the HWHM  remains finite even in the -CSI limit, e.g., at  = 0 and /2 of  XY , where we ought to have 1/  −→ 0 as  −→ 0. We believe that this is an artefact due to correlation cutoffs and other approximations used in PFFRG.To correct for it, we define the inverse correlation length 1/  by subtracting a constant background from HWHM  , with 1/ 0 equal to the HWHM  of -CSI.(D2) can also be understood as the statement that we are only interested in the change of HWHM  relative to that of -CSI as we move away from the Ising limit.We remark that the truncation range of  = 6 nearest neighbor bonds in our PFFRG calculations is greater than  0 , so subtracting the inverse of the former is not enough to yield a zero 1/  for -CSI.Other approximations in PFFRG, e.g.neglecting higher order vertices, might contribute to the background intensity as well.
The inverse correlation length calculated by (D2) is plotted as a function of  for the models (D1a) and (D1b) in Figs.11c  and 11d.These data suggest that we can interpolate smoothly from -QSI to -QSI or -QSI via an XY-type model or through the Heisenberg point, while keeping the density of monopoles associated with the dominant interaction small (less than one monopole every four cubic unit cells).Indeed, tracking the lowest of 1/  , we see that it only grows to approximately 0.1 × 2/ (at the Heisenberg point), which corresponds to a correlation length   ≈ 1.59 three to four times greater than the tetrahedral center-to-center distance √ 3/4 ≈ 0.43.The picture that emerges from this analysis is consistent with the scenario where the quantum spin liquids near the Ising limit, at the XX point, and at the Heisenberg point [83] are continuously connected in an extended -flux QSI phase.

FIG. 1 .
FIG. 1. Phase diagram of the  = 1/2 nearest neighbor XYZ model (2) on the pyrochlore lattice, in the parameter space where   = 1 and −1 ≤   ,   ≤ 1, obtained by the pseudofermion functional renormalization group.The labels 0-QSI, -QSI, Z(X)-AIAO represent the 0-flux quantum spin ice, the -flux quantum spin ice, and the all-in-all-out magnetic order along the local () axes, respectively.The empty squares indicate the parametrizations of Ce 2 Zr 2 O 7 proposed by Refs.[16,36] at the level of nearest neighbor interactions.The black dots demarcate the symmetry-preserving and symmetrybreaking regions as signalled by the breakdown of smooth renormalization group flow, with the error bars mainly reflecting the resolution of the grid for the calculations.The lines (without dots) separating the two spin ices and the two magnetic orders are obtained by further examining the diagonal components of the static susceptibility.See Sec.IV A for more details.

(
FIG. 6.The  component of the static susceptibility, and the spin-flip and non-spin-flip neutron scattering structure factors resulting from it, at (  ,   ) equal to (a-c) (1, 0.17) and (g-i) (0.63, 0.25), calculated using gauge mean field theory.The  component of the equal-time spin structure factor, and the spin-flip and non-spin-flip neutron scattering structure factors resulting from it, at (  ,   ) equal to (d-f) (1, 0.17) and (j-l) (0.63, 0.25), calculated using gauge mean field theory.Both parameters stabilize the -flux QSI.

FIG. 7 .
FIG. 7. Diagrammatic representations of the flows of (a) self-energy and (b) two-particle vertex, see (A5a) and (A5b).The slashed propagator in (a) represents the single-scale propagator  Λ () defined in (A2), while the pair of slashed propagators in (b) represents the sum of the two cases where one of the propagators is replaced by  Λ kat () defined in (A3).Site indices are conserved along solid lines.This figure is adapted from Ref. [57].

HWHMFIG. 11
FIG. 11.(a) We calculate the half width at half maximum (HWHM) for each diagonal component of the static susceptibility along the cut k = (0, 0, ) through the pinch point (see the inset where the cut is indicated by the white line).The example shown here is    (k) of the XY model (D1a) at / = 0.25, normalized by the maximum at    (k = 2X).(b) HWHM  extracted from   (k) for  = , ,  as a function of  for the XY model (D1a).(c) The inverse correlation length 1/  , calculated from HWHM  by subtracting a constant background (D2), as a function of  for the XY model (D1a).(d) The inverse correlation length 1/  as a function of  for the XXZ model (D1b).The arrows in (b-d) indicate divergences of HWHM , and 1/ , in the Ising limit (  ,   ,   ) = (1, 0, 0), where  , (k) are completely flat.
. Previous analyses on Ce 2 Zr 2 O 7 propose that   is dominant and antiferromagnetic, including the possibility of   ≈   [16, 18, 36].Here we shall set   = 1 as the unit for energy, and study the XYZ model in the parameter regime where |  |, |  | ≤   .Further considerations about the general parameter regime   ,   ,   ≥ 0 are given in Appendix D.