Ground state phase diagram of dipolar-octupolar pyrochlores

The dipolar-octupolar pyrochlore oxides R$_2$M$_2$O$_7$ (R=Ce, Sm, Nd) represent an important opportunity in the search for three dimensional Quantum Spin Liquid (QSL) ground states. Their low energy physics is governed by an alluringly simple XYZ Hamiltonian, enabling theoretical description with only a small number of free parameters. Meanwhile, recent experiments on Ce pyrochlores strongly suggest QSL physics. Motivated by this, we present here a complete analysis of the ground state phase diagram of dipolar-octupolar pyrochlores. Combining cluster mean field theory, variational arguments and exact diagonalization we discover multiple $U(1)$ QSL phases which together occupy a large fraction of the parameter space. These results give a comprehensive picture of the ground state physics of an important class of QSL candidates and support the possibility of a $U(1)$ QSL ground state in Ce$_2$Zr$_2$O$_7$ and Ce$_2$Sn$_2$O$_7$.

What has yet to be achieved is the combination of a material with an experimentally robust QSL state, with theoretical understanding of the microscopic interactions which give rise to that state and what kind of spin liquid they produce. Some materials studied as potential QSLs actually order at low temperature [17][18][19] , and others are complicated by chemical or structural disorder [20][21][22][23] . Meanwhile, the relevant theoretical models are often complicated, possessing many free parameters [24][25][26] .
In DO pyrochlores, the magnetic rare earth ions form a corner-sharing tetrahedral structure [ Fig. 1 (inset)]. There are strong crystal electric fields (CEFs) acting on each magnetic site, resulting in a Kramers doublet at the bottom of the CEF spectrum, separated from higher states by a large gap ∆ CEF ∼ 100K 30,37,40 . With the scale of exchange interactions being ∼ 1K 28,39 , this motivates a description of the system in terms of pseudospin-1/2 operators τ x i , τ y i , τ z i . The thing which sets DO pyrochlores apart from other pyrochlore oxides is the transformation properties of these operators under time-reversal and lattice symmetries 45,46 . τ x i and τ z i both transform like the component of a magnetic dipole ori-ented along the site's C 3 symmetry axis, while τ y i transforms like a component of the magnetic octupole tensor.
Assuming nearest-neighbor interactions, symmetry constrains the Hamiltonian to take the form 45 : The final term in Eq. (1) can be removed by a suitably chosen global transformation τ α →τα 45,48 , reducing the problem to an XYZ Hamiltonian: An understanding of dipolar-octupolar pyrochlores and their potential to realize QSL ground states requires understanding of the ground state phase diagram of Eq. (2). Certain limits of the parameter space of Eq. (2) have been well studied, namely: the perturbative limit where one exchange parameter dominates the other two 4,47,49 , the XXZ limit where two of the three exchange parameters are equal 5,50-54 and the region of parameter space without a sign problem for Quantum Monte Carlo (QMC) 5,45,[53][54][55] . However, there is no reason to expect materials of interest to fall into one of these limits, so a global phase diagram is needed.
In this Article, we calculate the ground state phase diagram of Eq. (2), by combining Cluster Mean Field Theory (CMFT), a variational extension to CMFT (CVAR) 52 and Exact Diagonalization (ED). Where the results can be compared with available QMC results 55 , they agree well. The final result for the phase diagram is shown in Fig. 1, with the parameter space expressed in terms of an overall scaleJ which can be divided out and two angles φ, ψ: Jx =J cos(φ) sin(ψ),Jỹ =J sin(φ) sin(ψ), We find four U (1) spin liquid phases, occupying a large combined portion of the parameter space, competing with an antiferromagnetic "all in/all out" (AIAO) phase Jx =J cos( ) sin( ),Jỹ =J sin( ) sin( ),Jz =J cos( ) < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 g p M 9 w A o R g V q 8 O t 6 f V E 1 r C a N e s A = " > A A A C g H i c b Z F b S 8 M w G I b T e p r 1 V P X S m + A Q N p i z n Y I i C E N v Z F c K 7 g D r G G m W b W F p W p J U r K W / w / / l n T 9 G M N t 6 4 Q 4 f B B 7 e 7 3 3 5 k i 9 + x K h U j v N j m B u b W 9 s 7 h V 1 r b / / g 8 M g + P m n J M B a Y N H H I Q t H x k S S M c t J U V D H S i Q R B g c 9 I 2 5 8 8 T f v t d y I k D f m b S i L S C 9 C I 0 y H F S G m p b 3 9 5 i r I B S R t Z P 8 3 x I 8 s e U s 9 H A j Y y 6 O F Q l r x o T M v Q k 5 R r l L R c 8 a z V W L I Q m 3 u X Y 3 A 1 9 r l m m v b 2 7 a J T d W Y F V 8 H N o Q j y e u n b 3 9 4 g x H F A u M I M S d l 1 n U j 1 U i Q U x Y x k l h d L E i E 8 Q S P S 1 c h R Q G Q v n S 0 w g x d a G c B h K P T h C s 7 U / 4 k U B V I m g a + d A V J j u d y b i u t 6 3 V g N 7 3 o p 5 V G s C M f z Q c O Y Q R X C 6 W / A A R U E K 5 Z o Q F h Q f V e I x 0 g g r P S f W X o J 7 v K T V 6 F V q 7 r X 1 d r r T b H + m K + j A M 7 A O S g B F 9 y C O n g G L 6 A J M P g 1 i k b F u D R N s 2 R e m e 7 c a h p 5 5 h Q s l H n / B 0 E I w r Q = < / l a t e x i t > AIAO Octupolar Order (2)] on the pyrochlore lattice (inset), describing dipolar-octupolar pyrochlores. The three exchange parametersJα are represented in terms of an overall scaleJ and two angular variables φ, ψ [Eq. (3)]. The phase diagram features "all in/all out" (AIAO) and octupolar ordered phases, and four distinct U(1) QSLs. These four QSLs are distinguished by whether the emergent electric field of the low energy gauge theory transforms like a magnetic dipole or octupole, and by the flux penetrating elementary plaquettes in the ground state (0 or π). The phase diagram is obtained by combining Cluster Mean Field Theory (CMFT), a cluster variational (CVAR) calculation and Exact Diagonalization (ED) as described in the text. The two regions bounded by black dashed lines correspond to the subset of parameters in Eq. (5), from which the entire phase diagram can be generated using unitary transformations. and octupolar order. The four U (1) QSLs all host gapless photons and gapped fractionalized charges, and are thus realizations of emergent electromagnetism 4,6,47,49 . They are labelled dipolar/octupolar-U (1) 0/π with the dipolar/octupolar label referring to whether the emergent electric field transforms like a magnetic dipole or octupole 56,57 , and the 0/π subscript referring to the U (1) flux penetrating elementary plaquettes in the ground state.
The remainder of this Article is devoted to explaining the calculations leading to Fig. 1, before finishing with a brief discussion of the outlook for experiments.
The Article is structured as follows: • In Section II we describe some simple dualities which allow the whole phase diagram to be generated from calculations covering only a subregion of parameter space.
• In Section III we calculate the ground state phase diagram using CMFT, augmented with the CVAR approach.
• In Section IV we show ED calculations on a 16site cluster, and use these as an alternative route to calculate the ground state phase diagram.
• The construction of the complete phase diagram [ Fig. 1], from the combination of the calculations in the preceding sections, is then described in Section V.
• Section VI gives a summary of the results and an outlook for future work on dipolar-octupolar pyrochlores.

II. DUALITIES OF THE MODEL AND REDUCED PARAMETER SPACE
In calculating the phase diagram it is useful to note that Eq.
where U γ,v represents a global rotation by an angle γ around axis v of pseudospin space (which is not the same as a rotation in the physical crystal space). Making use of these dualities means that we do not actually need to (2), (6)] within the region of parameter space given by (7) withJz > 0. CMFT calculations give two regimes for the optimal configuration of the auxiliary fields hi: an ordered region where the hi point uniformly along the y axis of pseudospin space (red) and a disordered region with a large degeneracy of CMFT solutions where hi = σihzi , with signs σi summing to zero on every tetrahedron (green). The CVAR calculation, which incorporates quantum tunnelling between CMFT solutions, breaks the degenerate region into two, based on the sign of the effective tunnelling matrix element g ef f . Positive (negative) values of g ef f lead ultimately to a π-flux (0-flux) U(1) QSL ground state. study the full parameter space ofJx,Jỹ,Jz, it is enough to consider a subset of parameters |Jz| > |Jx|, |Jỹ|,Jx >Jỹ (5) from which we can then generate the rest of the phase diagram by applying the transformations from Eq. (4) to our results. The parameter space dsecribed by (5) is delineated by the black dashed lines in Fig. 1. TakingJz to be the strongest exchange parameter as in (5), ifJz < 0 it is clear that the ground state will simply order ferromagnetically with respect to thez-axis of pseudospin space. In terms of the physical magnetic moments this implies AIAO order. The more challenging problem is to discover what happens whenJz > 0.
To study this case we rewrite the Hamiltonian in terms of spin ladder operatorsτ ± i : Jx +Jỹ andJ ±± = 1 4 Jx −Jỹ . The subregion of parameter space given by (5) then becomes:

A. CMFT Calculation
To begin, we consider the phase diagram using a tetrahedral CMFT, as employed for the XXZ limit (J ±± = 0) in Ref. 52. A summary of the calculation is given here, with a detailed description found in Appendix A.
To construct the CMFT we use the fact that the pyrochlore lattice can be divided into two sets of tetrahedra 'A' and 'B', with all neighbors of an 'A' tetrahedron being 'B' tetrahedra and vice versa. We then seek to optimize a product wave function over all 'A' tetrahedra: The wave function |φ t on each tetrahedron t is defined to be the ground state of a single tetrahedron Hamiltonian H t contains the original exchange terms acting on the bonds of t as well as auxiliary fields h i on each site The auxiliary fields h i then serve as variational parameters for optimizing |ψ CMFT , and a CMFT wave function can be indexed by a configuration of h i on the lattice. There are two regimes for the optimal configuration of h i in CMFT as shown in Fig. 2. For sufficiently large, positive, values ofJ ± orJ ±± the optimal solutions have h i ordered ferromagnetically along the y-axis of pseudospin space. This implies τỹ = 0, and therefore octupolar order sinceτỹ transforms like a magnetic octupole 45 .
In the remainder of the phase diagram there is a large, ice-like, degeneracy of disordered CMFT solutions, with h i = σ i hz i where h is a fixed, uniform, magnitude and σ i = ±1, subject to the constraint that σ i sum to zero on every tetrahedron.

B. Cluster variational (CVAR) calculation
To resolve the CMFT degeneracy in the disordered regime, we follow the cluster variational (CVAR) method 52 . The calculation is described briefly here with further details given in Appendix B. . The gap collapses rapidly upon crossing the white line, supporting the conclusion that this line corresponds to a transition to long range order breaking π-rotation symmetry around thez-axis in the thermodynamic limit.
Labelling CMFT ground states according to their configuration of signs {σ} we write down a generalized superposition of CMFT solutions where a {σ} are unknown coefficients. We then seek to optimize the new variational energy Eq. (12) can be expanded in terms of the overlap between distinct CMFT wavefunctions, in a similar spirit to the derivation of dimer models from an expansion in the overlap between singlet coverings of a lattice 59 . This generates an effective Hamiltonian in the space of CMFT solutions, where the leading term is a six-site ring exchange which flips the values of σ i on hexagonal plaquettes where σ alternates in sign around the plaquette, with matrix element g ef f . This Hamiltonian has already been studied using Quantum Monte Carlo 47,49 . It can have two different QSL ground states depending on the sign of g ef f . Both are U(1) QSLs with gapped, bosonic, charges and gapless photons. The two ground states are distinguished by the U(1) flux threading elementary plaquettes in the ground state. This background flux vanishes for g ef f < 0 (U(1) 0 ) but is equal to π on every plaquette for g ef f > 0 (U(1) π ). The value of g ef f can be extracted from the CMFT calculation for all values of exchange parameters (see Appendix B), and by this means the degenerate region within CMFT can be divided into two ground state QSL phases (U(1) 0 and U(1) π ) depending on the sign of g ef f . The boundary between regions with different signs of g ef f is shown in Fig. 2. This constitutes our estimate of the boundary between 0-flux and π-flux QSLs.

IV. EXACT DIAGONALIZATION
We now turn to ED calculations on a 16-site cubic cluster with periodic boundaries, to obtain alternative estimates of the phase boundaries.
A. Boundary of octupolar ordered phase Fig. 3(a) shows the second derivative of the ground state energy on this cluster with respect toJ ± , at various fixed values ofJ ±± /Jz. This second derivative exhibits a peak asJ ± is swept, indicating a qualitative change in the ground state 58 . Fig. 3(b), shows the position of these peaks as a function ofJ ± /Jz andJ ±± /Jz, laid over a color plot of the gap to excitations with odd totalτz. The parity p = (−1) iτz i is conserved by H, with the ground state always having p = 1. The line of peaks in the second derivative of the ground state energy coincides with a rapid decrease of the gap to p = −1 excitations. This suggests the formation of a twofold degenerate ground state in the thermodynamic limit, breaking π rotation symmetry around thez axis, consistent with the octupolar order identified in CMFT. We thus interpret the peaks in the second derivative of the ground state energy as indicative of a transition to octupolar order.

B. Boundary between QSLs
It is not easy to cleanly distinguish between the two QSL phase, QSL 0 and QSL π using ED on a small cluster. However, some insight into how to identify the phase boundary can be gained by considering how this transition occurs in the perturbative limit and in the CVAR approach.
From the perspective of both perturbation theory and CVAR, the transition from QSL 0 to QSL π occurs when the leading tunnelling matrix element, g ef f , between ice-like states changes sign. At the point where g ef f vanishes, tunnelling is restricted to higher order processes and will therefore be suppressed, leading to a near restoration of the degeneracy of ice-like states.
Returning to ED, this suggests that the transition from QSL 0 to QSL π will be accompanied by a simultaneous collapse of many excited states, in the sector with even totalτz, to near zero energy. Such a collapse is indeed observed in the ED data, as shown in Fig. 4(a). The position of this collective minimum in the gaps within the even sector constitutes the ED estimate of the phase boundary between the two QSLs. The phase boundary thus obtained is compared with that from CVAR in Fig.  4(b), with the two estimates agreeing closely.

C. Combining information from CMFT/CVAR and ED
Combining the information from CMFT/CVAR and ED gives the phase diagram shown in Fig. 5.
ForJ ± < 0 the CMFT and ED estimates of the octupolar phase boundary agree closely. ForJ ± > 0 the ED estimates a larger region of octupolar order (and hence a smaller QSL region) than does the CMFT approach.
ForJ ± > 0 the model has no sign problem from the perspective of QMC, and in this regime we can compare with previous QMC studies. Several previous QMC studies of the caseJ ±± = 0 have observed the transition from QSL 0 to the ordered phase asJ ± is increased 5,53-55 . A recent QMC study by Huang et al 55 has studied the behavior of this phase boundary as a function ofJ ±± . Comparison with these results can be used to adjudicate between ED and CVAR where they disagree. The ED calculation gives closer agreement with the QMC results from [55] than CMFT/CVAR does, and therefore we will take the ED calculation as our estimate of the boundary of the octupolar phase.
The estimates of the boundary between the two QSL phases agree closely between CVAR and ED, as shown in Fig. 4 (b). There is, however, some difference between the two estimates at larger negative values ofJ ± . For the purpose of Fig. 5 we use the boundary from CVAR because it gives a more direct prediction of the transition between the two states in the thermodynamic limit, as opposed to the more indirect inference from the behavior of gaps in ED.

V. CONSTRUCTION OF COMPLETE PHASE DIAGRAM
The phase diagram in Fig. 5 can then be extended to the full parameter space using the duality relations [Eq. (4)].
In doing this, we must take into account how the duality transformations act on the ground states. For exam-  (7) withJz > 0., obtained from combining CMFT, CVAR and ED calculations. This region of parameter space corresponds to the lower region bounded by dashed lines in Fig. 1. The full phase diagram in Fig. 1 can be generated by applying the dualities described in Sec II to this region and to the upper region bounded by dashed lines in Fig. 1 which has all-in-all-out order throughout.
ple, the octupolar ordered phase with τỹ i = 0 becomes an AIAO phase when acted on by a transformation which swaps theỹ-axis with thex-axis orz-axis. On the other hand, transformations which swap only thex-axis and z-axis, don't change the classification of the ground state phase becauseτx i andτz i transform equivalently under point-group and time reversal symmetries.
Similar considerations allow us to distinguish four different kinds of U(1) QSL, generated from the two in the phase diagram of Fig. 5. In the 0-flux and π-flux QSLs in Fig. 5 the emergent electric field of the QSL E i ∼τz i

4,50
and therefore transforms like a magnetic dipole. If we act a transformation that swaps thez-axis andỹ-axes, then E i ∼τỹ i and transforms like an octupole. We should therefore not only distinguish U(1) QSLs by the flux but by the dipolar or octupolar character of the emergent electric field, giving four distinct QSLs on the complete phase diagram 56,57 .
Applying these arguments to the parameter space covered in Fig. 5, and to the case ofJz < 0, allows us to generate the full phase diagram, shown using spherical coordinates [Eq. (3)] in Fig. 1.

VI. SUMMARY AND OUTLOOK
We have thus established a phase diagram for the generic, symmetry allowed, nearest neighbour exchange Hamiltonian describing dipolar-octupolar (DO) pyrochlores R 2 M 2 O 7 (R=Ce, Sm, Nd). The picture we arrive at is an encouraging one for the realization of QSL states. There are four distinct U(1) QSLs on the phase diagram of the generic nearest neighbor model, and between them they occupy ∼ 19% of the available parameter space.
Amongst materials, Ce 2 Zr 2 O 7 29,30 , Ce 2 Sn 2 O 7 27,28 and Sm 2 Zr 2 O 7 35 stand out as lacking low temperature order. The Ce pyrochlores in particular seem promising with recent neutron scattering results on Ce 2 Zr 2 O 7 bearing similarity to predictions for emergent photons 30 . Low energy correlations in Ce 2 Sn 2 O 7 seem to be dominantly octupolar in nature 28 , which would be consistent with either of the two octupolar spin liquids on the phase diagram [ Fig. 1].
It will be important to establish estimates of the exchange parameters of Ce 2 Zr 2 O 7 and Ce 2 Sn 2 O 7 , combining information from inelastic neutron scattering with fits to thermodynamic data. Mean field calculations in [28] give an initial estimate for Ce 2 Sn 2 O 7 of J y = 0.48K, J z = 0.03K, while setting J x and J xz to zero, in the basis of Eq. (1). This would place Ce 2 Sn 2 O 7 in the Octupolar-U (1) π region of the phase diagram. It would be useful to refine this estimate with all parameters allowed to be finite, and using calculations beyond mean field theory.
If refined parameterisations place Ce 2 Zr 2 O 7 and Ce 2 Sn 2 O 7 within the QSL regimes of Fig. 1, then this will be a strong indication that they are indeed U(1) QSLs, and the parameterized model will provide a platform for further theoretical study. Understanding the effects of disorder of the crystal structure is also likely to be crucial, particularly in regard to the possible substitution of magnetic Ce 3+ with non-magnetic Ce 4+30 .
For those DO pyrochlores that are known to possess magnetic order at low temperature, the spin liquid phases may also manifest at finite temperature, as suggested recently in Nd 2 Zr 2 O 7 61 . In such cases it may even be possible to tune into the T = 0 QSL phase using chemical or physical pressure, giving another avenue to realize these exotic states of matter. The CMFT proceeds by optimizing variational wavefunctions of the form: where the product is over all 'A' tetrahedra [ Fig. 6], {h i } is the configuration of auxiliary fields defined on each site and the single tetrahedron wavefunctions |φ t depend only on the fields on sites belonging to tetrahedron t. The auxiliary fields h i are variational parameters for optimizing the CMFT energy The wave functions |φ t ({h i∈t }) are taken to be eigenstates of a single tetrahedron Hamiltonian H t E CMFT is then: where the final term in Eq. (A5) sums over bonds belonging to 'B' tetrahedra and accounts for the interactions on those tetrahedra. There is a large region of the phase diagram [ Fig. 2 of main text] in which the optimal solutions for h i take the form: Correspondingly, the expectation values of the spin components are: with h and s being uniform across the system, and fixed by the energy optimization for a given parameter set. With this form for the auxiliary fields, the mean field energy [Eq. (A5)] becomes Any arrangement of signs σ i such that gives rise to the same value of 0,t , as can be inferred from the symmetries of the original Hamiltonian. The remaining terms in Eq. (A8) are also the same for all configurations obeying Eq. (A9). Thus we have a large degeneracy of mean field solutions in this regime. Each arrangement of signs σ i obeying Eq. (A9) defines a CMFT wavefunction [via Eqs. (A1), (A4) and (A6)] which we will denote with |ψ CM F T ({σ}) . Explicitly, the form of single tetrahedron wave functions |φ t ({σ i∈t ) (denoted simply as |σ 0 σ 1 σ 2 σ 3 ) relates to the configuration of signs on t in the following way, written in the basis diagonalizingτz i : The parameters µ, ν and ρ can always be chosen to be real. This choice, combined with the choice to define the first term on the right hand side of each line of (A10) to be positive, removes any phase ambiguity in the CMFT wavefunctions. µ, ν and ρ vary as a function of the exchange parametersJα and are plotted in Fig. 7.

Appendix B: Details of CVAR calculation
The goal of the CVAR calculation is to resolve the degeneracy of the CMFT solutions by considering a new trial wavefunction which is a superposition of the CMFT solutions: where a {σ} are, a priori unknown, complex, coefficients. We then seek to optimize the variational energy where X is a matrix containing the Hamiltonian matrix elements between different CMFT wavefunctions and O contains the overlaps (the CMFT wavefunctions are not generally orthogonal to one another) It is then useful to define a new matrix X with vanishing diagonal elements: such that We then relate the vector of coefficients a, to a new normalized vector b via: The variational energy is then The optimal superposition of CMFT solitions is then given by the ground state of H eff and Eq. (B7). We then expand H eff in terms of two overlap parameters o 2 and o 4 , which can be defined from the wavefunctions in Eqs. (A10): These two quantities are treated as small parameters for the purposes of the expansion and indeed they are small through most of the relevant parameter space, as shown in Fig. 8. To expand Eq. (B10) we note that all the diagonal elements of O are unity, and the leading off diagonal elements ∼ (o 2 ) 3 (coming from the process illustrated in Fig. 9) so we can write The first two terms of the expansion of H eff are then: The leading elements in X are ∼ (o 2 ) 2 (again, from the process in Fig. 9) and the leading elements in O · X are ∼ (o 2 ) 5 and so we henceforth drop the second term. We then need to evaluate the leading matrix elements in X which connect configurations of σ i which differ on a single hexagonal plaquette as shown in Fig. 9. The matrix element to a flip a hexagon is g ef f . The sign of g ef f determines whether the ground state should be a 0 or π flux QSL, with as may be inferred from prior quantum Monte Carlo studies of the six-site ring exchange Hamiltonian 49 and from a unitary transformation which relates the sign-problem free case (g ef f < 0) to the frustrated case (g ef f > 0) 4 .
Quite generally the matrix element of X between two CMFT wavefunctions can be written as where the first sum is over 'A' tetrahedra and the second is over bonds belonging to 'B' tetrahedra. H t is the original exchange Hamiltonian on tetrahedron t (distinct from H t in Eq. (A3)) and −2J ± µ(µ + 2(ν + 1 − µ 2 − ν 2 − ρ 2 )). (B18) The contribution of any 'A' tetrahedron which does not change configurations between {σ } and {σ} to Eq. (B17) vanishes. Similarly, the contribution of any 'B' bond connecting two unchanged 'A' tetrahedra vanishes.
Secondly, there are contributions from 'B' bonds connecting to the interior of the flipped hexagon (highlighted in blue in Fig. 10): Finally, there are contributions from 'B' bonds connecting to the exterior of the 'A' tetrahedra on the flipped hexagon (highlighted in green in Fig. 10): Summing these contributions and accounting for the fact that σ i must alternate around the hexagon and must obey Eq. (A9) everywhere, we arrive at the matrix element: From Eq. (B24) and the numerically determined CMFT wavefunctions we can calculate g ef f and thus predict the ground state in the degenerate region of CMFT from the sign of g ef f . The behavior of g ef f and sign(g ef f ) over the relevant region of parameter space is shown in Fig. 11.