title: Hierarchy of energy scales and field-tunable order by disorder in dipolar-octupolar pyrochlores
authors: Benedikt Placke Roderich Moessner Owen Benton
text: 
  
  
  
  
  
  
  Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Straße 38, 01187 Dresden, Germany
  
  
  
  
  
  
  
  Dipolar-octupolar pyrochlore magnets in a strong external magnet field applied in the [110] direction are known to form a `chain' state, with subextensive degeneracy. Magnetic moments are correlated along one-dimensional chains carrying effective Ising degrees of freedom which are noninteracting on the mean-field level.
  Here, we investigate this phenomenon in detail, including the effects of quantum fluctuations.
  We identify two distinct types of chain phases, both featuring distinct subextensive, classical ground state degeneracy. 
  Focussing on one of the two kinds, we discuss lifting of the classical degeneracy by quantum fluctuations. We map out the ground-state phase diagram as a function of the exchange couplings, using linear spin wave theory and real-space perturbation theory. 
  We find a hierarchy of energy scales in the ground state selection, with the effective dimensionality of the system varying in an intricate way as the hierarchy is descended. We derive an effective two-dimensional anisotropic triangular lattice Ising model with only three free parameters which accounts for the observed behavior. 
  Connecting our results to experiment, they are consistent with the observation of a disordered chain state in Nd_2Zr_2O_7. We also show that the presence of two distinct types of chain phases has consequences for the field-induced breakdown of the apparent U(1) octupolar quantum liquid phase recently observed in Ce_2Sn_2O_7.
  
  
  
  
  
  Hierarchy of energy scales and field-tunable order by disorder in dipolar-octupolar pyrochlores
      Owen Benton
      December 30, 2023
  ===============================================================================================
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Frustrated magnets consisting of corner-sharing fully connected units often feature nontrivial ground state degeneracy <cit.>.
  This is captured by local constraints on the configuration of each unit, which nevertheless leave the system with a large number of unconstrained degrees of freedom. In the case of the pyrochlore lattice, which consists of corner-sharing tetrahedra, this often takes the form of an ice-like two-in-two-out rule, leading to the class of materials known as “spin-ice” <cit.>.
  
  One natural question to ask is what is the fate of this classical degeneracy in the presence of quantum fluctuations. Since the degeneracy is “accidental” in the sense that the degenerate states are not related by a symmetry of the Hamiltonian, it is expected to be lifted by fluctuations, a process which is generally called order by disorder (OBD) <cit.>. Probably the most prominent experimental example of quantum OBD is the effective spin-12 pyrochlore Er_2Ti_2O_7 <cit.>. 
  In Er_2Ti_2O_7 the OBD mechanism lifts the degeneracy of a manifold of antiferromagnetic states with a single U(1) degree of freedom. There are competing contributions to the ground state selection <cit.> which cannot be easily tuned, making them difficult to disentangle from one another.
  
  In this article, we consider a different, experimentally motivated, case of OBD on the pyrochlore lattice: the case of dipolar-octupolar pyrochlores in a [110] magnetic field.
  Here, in contrast to Er_2Ti_2O_7, the mean field degeneracy grows sub-extensively.
  Both the energy scale of ground state selection, and the selected ground state can be tuned by an external field, as also seen in some theoretical models on other lattices <cit.>.
  
  Dipolar-octupolar pyrochlores are a family of rare-earth oxides R_2M_2O_7 with R=Nd, Ce, Sm <cit.>.
  They are well described by a nearest-neighbor XYZ model with a peculiar coupling to an external magnetic field <cit.>. The dipolar-octupolar crystal field doublet couples only to the projection of the external field onto the local [111] easy axes, which are not collinear. A uniform field applied in the crystalline direction [110] thus separates the lattice into two sets of one-dimensional chains (see <ref>), with a frustrated interaction that cancels out exactly on the mean-field level. 
  
  Here, we investigate this phenomenon through calculations of the classical phase diagram and by showing how quantum fluctuations lift the degeneracy of the chain states and restore correlations between chains, albeit at extremely small energy scales.
  
  Studying the classical phase diagram as a function of exchange parameters and field, we find four distinct chain phases, which can be divided into two different types. First, there are three Chain_λ (λ = x, y, z) phases, in which the chains along the field direction (α chains) are polarized, while those perpendicular to it (β chains) each retain one Ising-like degree of freedom [see <ref>]. This yields a subextensive ground-state entropy scaling as the square of the linear system size L^2 when N∼ L^3. 
  
  There is then one further type of chain phase, Chain_Y^* which appears due to the special status of the pseudospin y-axis in the XYZ model for dipolar octupolar pyrochlores. The operator S^y_i transforms like a magnetic octupole <cit.> and therefore does not couple linearly to an external field. This gives rise to the Chain_Y^* phase at low fields if the octupolar exchange coupling J_y is dominant.
  At low fields, the moments on the α chains  have a finite octupolar component, which results in an additional independent Ising degree of freedom per α chain.
  The existence of this additional chain phase for strong J_y implies the possibility of a field-induced confinement transition out of the octupolar quantum liquid state recently proposed for Ce_2Sn_2O_7<cit.>.
  
  Turning to the effect of quantum fluctuations, we show that for
  the Chain_λ phases, fluctuations via the polarized α chains mediate an effective interaction between the unpolarized β chains and thus lift the subextensive classical degeneracy. We compute the leading effective chain interactions and map out the ground-state phase diagram in the semiclassical limit S→∞. 
  To this end, we use real space perturbation theory (RSPT) <cit.> to derive a two-dimensional model describing the effective interaction between the chains, which takes the form of a nearest-neighbor anisotropic triangular lattice Ising model.
  To corroborate our results, we also compute the zero point energy of different classical ground state spin configurations numerically using linear spin wave theory (LSWT) and find that the full splitting of the classically degenerate states can be captured by the effective model remarkably well. The two methods agree also on a quantitative level, with positions of phase boundaries agreeing even for moderately strong transverse couplings.
  
  Two kinds of ground state are selected by OBD, depending on the exchange parameters and external field. First, a completely uniform state, where all β chains align ferromagnetically and second, zigzag configurations where β chains align antiferromagnetically in the [110] direction.
  The selection is governed by two competing effective nearest-neighbor chain interactions generated by correlated fluctuations around hexagons. The two couplings scale differently as a function of the external magnetic field, making it possible, in some parts of the phase diagram, to tune between the ferromagnetic and zigzag ground states.
  The left-over degeneracy in the zigzag phase is lifted on a minuscule energy scale that however can be resolved in our LSWT calculation. In terms of the effective model, the splitting can be fully understood by including an effective next-nearest-neighbor coupling, which appears at higher order in RSPT.
  
  Discussing the effect of temperature we compute the classical low-temperature expansion around different ground states. Remarkably, we find that, in contrast to its quantum counterpart, at leading order classical OBD is completely absent in all Chain_λ phases.
  
  Guruciaga et al.<cit.> have considered a related problem of thermal OBD in a classical pyrochlore Ising model in a [110] field.
  The case considered here differs from [Guruciaga2016] in that we consider quantum OBD with a Hamiltonian directly motivated by experiments on dipolar-octupolar pyrochlores<cit.>, and in that we find multiple different ground states as a function of field and exchange parameters.
  
  The remainder of this paper is organized as follows. In <ref>, we review dipolar-octupolar pyrochlores and their classical ground-state phase diagram in a [110] field.
  This serves as a starting point for the OBD calculation presented in the following. In <ref>, we derive a two-dimensional triangular lattice Ising model describing the lifting of the classical degeneracy in leading-order real space perturbation theory (RSPT). In <ref>, we compute the ground-state phase diagram using linear spin wave theory (LSWT), comparing the
  results with those from RSPT.
  In <ref> we discuss the tunability of the ground state selection as a function of the external field strength. We compare the results of the LSWT calculation with those from RSPT. Section <ref> explores the connection of our work to experiment. We conclude in <ref>.
  
  
  
  
  
  
  
  § DIPOLAR-OCTUPOLAR PYROCHLORES
  
  
  
  
   §.§ Exchange Hamiltonian
  
  
  When describing the low-energy properties of rare-earth magnets, we have to start with the the local crystal electric
  field (CEF). The CEF energy scales are typically orders of magnitude larger than the exchange energy between two ions <cit.>. To calculate low temperature, low energy properties, it is therefore usually sufficient to consider the lowest energy multiplet of the crystal field. For materials within our present interest, this takes the form of a doublet.
  
  On the pyrochlore lattice [<ref>] this doublet comes in three different flavors.
  Each flavor corresponds to an irreducible representation of the double group of the point group D_3d, together with time reversal invariance. In the case that the two states forming the lowest CEF are not related by spatial symmetry but only by time-reversal, it is called a dipolar-octupolar doublet <cit.>. This occurs in rare-earth pyrochlores based on Nd, Ce and Sm.
  
  Once projected into a dipolar-octupolar doublet, the only surviving component of the magnetic moment μ⃗_i operator is that oriented along the local [111] easy axis z_i, which connects the centers of the two tetrahedra sharing site i. In our convention, the z_i are constant on each of the four sublattices ℒ_ν (ν=0, 1, 2, 3) and point into the tetrahedron shown in figure <ref> (b).
  
  Defining pseudospin-12 operators S^α_i acting within the doublet, one can choose a basis such that
  
      μ⃗_i = μ_B g_z  z_i S_i^z.
  
  
  While S^x_i and S^z_i transform like the z_i component of a magnetic dipole moment, S^y_i transforms like a component of a magnetic octupole tensor <cit.>. Expressed in terms of components of the angular momentum operator 𝒥, the corresponding octupole moment is <cit.>
  
      S^y∝𝒥_y^3 - 𝒥_y 𝒥_x^2 - 𝒥_x 𝒥_y 𝒥_x - 𝒥_x^2 𝒥_y,
  
  which transforms trivially under all lattice symmetries (D_3d) but is odd under time reversal.
  
  
  The most general symmetry allowed nearest-neighbor exchange Hamiltonian for such a dipolar-octupolar pyrochlore magnet then only couples pseudospin components with the same transformation properties
  
  
      ℋ =   ℋ_ex + ℋ_field, 
      ℋ_ex =   ∑_ij J_xx S_i^x S_j^x + J_yy S_i^y S_j^y + J_zz S_i^z S_j^z 
         ∑_ij+ J_xz( S_i^x S_j^z + S_i^z S_j^x ),
      ℋ_field =     -μ_B g_z ∑_i H⃗· z_i S_i^z,
  
  
  where the first sum in the exchange Hamiltonian ℋ_ ex runs over all nearest-neighbor pairs ij in the pyrochlore lattice and the sum in the field Hamiltonian ℋ_ field runs over all sites i.
  
  We can remove the XZ exchange by a rotation of the pseudospin axes by an angle θ around the pseudospin y-axis
  
  
      S_i^x    →cosθ S_i^x - sinθ S_i^z, 
      
          S_i^z    →sinθ S_i^x + cosθ S_i^z, 
      tan(2θ)    = 2 J_xz^(0)/J_xx^(0) - J_zz^(0),
  
  
  which yields
  
  
      ℋ   = ℋ_ex + ℋ_field, 
      ℋ_ex   = ∑_ij J_x S_i^x S_j^x + J_y S_i^y S_j^y + J_z S_i^z S_j^z, 
      ℋ_field   =  - μ_B g_z ∑_i H⃗· z_i (sinθ  S_i^x + cosθ  S_i^z).
  
  
  
  To avoid too ornamented notation, we have not renamed the spin operators S_i^α. From here on, the S_i^α refer to those in <ref>.
  Altogether, dipolar-octupolar pyrochlores are described by a relatively simple XYZ exchange Hamiltonian, however with a peculiar coupling to the magnetic field.
  
  
  
   §.§ Classical phase diagram in a [110] field 
  
  
  
  
  Before turning to the case of a finite field, we quickly discuss the zero-field phase diagram. For zero field (H⃗=0) there are six possible classical ground state phases. 
  
  There are the three ice-like phases SI_λ  (λ=x,y,z), stabilized when the
  corresponding exchange coefficient is positive and sufficiently strong
  
      J_λ > max(-3J_λ', J_λ')   ∀λ' ≠λ.
  
  In these phases, all pseudospins align with the λ axis and are governed by a local “two up, two down” constraint with respect to this axis on every tetrahedron. This rule leaves an extensive number of degrees of freedom. If λ=x,z then the rule corresponds to
  a “two-in-two-out" rule on the magnetic moments; i.e. for each tetrahedron, two of the four magnetic moments on its corners point into and two point out of it. SI_x and SI_z thus describe the well known spin ice state, whereas SI_y describes a state analogous to spin ice but where the active degrees of freedom are octupolar moments [<ref>].
  
  There are then three ordered, “all-in-all-out” phases (AIAO_λ,  λ=x,y,z), occurring for sufficiently strong negative exchange coefficient
  
      J_λ < min(-J_λ'/3, J_λ')   ∀λ' ≠λ.
  
  These correspond to ferromagnetic order of the pseudospins along the λ axis. If λ=x,z then this will correspond to antiferromagnetic all-in-all-out order of the magnetic moments. The AIAO_y phase, by constrast, is a form of octupolar order.
  
  Turning to finite fields, for the rest of this article we set 
  
      H⃗ = √(3) h/2μ_B g_z(1
       1
       0).
  
  Substituting this into <ref> yields:
  
      ℋ_field =     - h ∑_i∈ℒ_0(sinθ  S_i^x + cosθ  S_i^z) 
         +h ∑_j∈ℒ_3(sinθ  S_j^x + cosθ  S_j^z),
  
  where ℒ_ν  (ν=0,1,2,3) denotes the four fcc sublattices, also indicated in <ref> (b). The field, when applied in this direction, only couples to sites on the sublattices 0 and 3. As shown in <ref> (a), sites on which the magnetic moment couples to the external field then form one-dimensional chains which are parallel to the field direction (α chains), whereas sites whose magnetic moment does not couple to the external field also form one-dimensional chains, which are perpendicular to the field direction (β chains). 
  
  
  The classical ground state phase diagram of dipolar-octupolar pyrochlores with H∥(1,1,0) is shown in  <ref>. There are now seven possible phases, six of which are directly connected to the six phases possible at zero field. The SI_Y phase evolves into two distinct chain phases as a function of field strength, separated by a second-order transition. 
  
  At any finite field the extensive degeneracy of the three spin ice phases SI_λ is broken down to a subextensive degeneracy, yielding effectively one-dimensional degrees of freedom.
  We first consider the dipolar cases, i.e. λ=x,z, shown in <ref> (b): while the magnetic moments on the α chains (sublattices 0 and 3) are polarized by the field, magnetic moments on the β chains do not couple to the field at all, but are restricted to two possible configurations per β chain by the two-in-two-out rule. Thus, each β chain carries an independent effective Ising degree of freedom. We denote this phase when resulting from dimensional reduction in a SI_λ phase by Chain_λ. In all of them, there is a finite magnetization in the direction of the field, stemming from the polarized α chains and each β chain carries a finite magnetization in the direction of the chain multiplied by its respective Ising degree of freedom.
  
  The case of the octupolar spin ice phase SI_Y is special since the ordered moment does not couple to the field. Because of this, the low-field phase which evolves out of the octupolar ice SI_Y is distinct from the other chain phases.
  If the octupolar exchange J_y dominates, a weak field partially polarizes the magnetic moments on the α chains, but leaves a finite octupolar component S_i^y. This component does not couple to the external field even on the α chains, however the “spin-ice” rule ∑_tet S^y_i = 0 enforces ordering of those moments along both α and β chains. 
  Hence, again the extensive degeneracy is broken down to a subextensive degeneracy, but with an Ising degree of freedom residing on each α and each β chain. We call this phase Chain_Y^*.
  
  
  
  The additional Ising degrees of freedom in the Chain_Y^* phase are related to the staggered octupolar moment on the α chains, which on a chain C is given by
  
      m_y^ Staggered = 1/2S(∑_i ∈ℒ_0 ∩ C  S_i^y -∑_j ∈ℒ_3 ∩ C S_j^y).
  
  As the field h is increased, or the octupolar exchange coupling J_y is lowered, this quantity vanishes in a second order transition [<ref> and <ref>] marking the onset of a chain phase Chain_Y which is equivalent to the dipolar chain phases Chain_X/Z in that there is exactly one Ising degree of freedom per β chain.
  
  The three AIAO_λ phases are stable in a small magnetic field. However, if the respective coupling J_λ is sufficiently weak (i.e. if  |J_λ| < J_λ'  for  some  λ' ≠λ), then the order is destroyed at a finite field, marking a transition to a Chain_λ' phase.
  Note that in a finite field, the AIAO_λ states will acquire a finite magnetization since the magnetic moments on the α chains will be partially polarized. However, even as h→∞, the phase retains an intensive two-fold degeneracy and a nonzero order parameter ∑_i S_i^λ≥ SN/2. These two AIAO states are related by a mirror symmetry in the case of the dipolar AIAO_X and AIAO_Z phases and by a combination of time-reversal and a mirror symmetry in the case of the octupolar AIAO_Y phase.
  
  
  
  
  
  § EFFECTIVE TRIANGULAR LATTICE ISING MODEL
  
  
  
  
  As established in the previous section, if none of the exchange constants J_λ are strongly negative, the system at moderate fields is in a chain phase [<ref>], where the classical ground state is subextensively degenerate with an Ising degree of freedom carried by each β chain [<ref>]. The interaction between those Ising degrees of freedom is frustrated and cancels out exactly on the mean-field level. However, since the resulting classical ground state degeneracy is accidental in the sense that the degenerate states are not related by any symmetry of the Hamiltonian [<ref>] one would expect it to be lifted by quantum fluctuations, a process called (quantum) order-by-disorder (OBD).
  
  
  
  
   §.§ Nearest-neighbor model
  
  
  In this section, we derive the leading order contribution to the effective interaction between the chains using real space perturbation theory (RSPT) <cit.>. The method divides the Hamiltonian into an unperturbed part, which has a particular classical state as its ground state, and a perturbation incorporating all transverse couplings. One then performs standard perturbation theory to obtain the corrections to the classical energy due to the transverse terms. Since the perturbation depends parametrically on the explicit classical state, the energy correction will depend on it too, possibly lifting classical degeneracies.
  
  We restrict the discussion to the “simple” chain phases Chain_λ where the effective interactions take the form of a nearest-neighbor anisotropic triangular-lattice Ising model.
  
  We can parameterize any classical ground state {S⃗_i^(0)} by L^2 Ising variables {η_i}, η_i ∈{+1,-1}, where L is the linear system size and N=4L^3 the number of sites in the lattice. We then want to calculate the lowest-order energy correction that distinguishes different ground states, that is different configurations of {η_i}. As we discuss in detail in app:rspt,
  such corrections arise from linked clusters of nontrivial topology on the pyrochlore lattice the smallest of which are hexagons [Fig. <ref> (a)]. The leading correction distinguishing configurations of {η_i} thus occurs at sixth order in RSPT.
  
  Any hexagon on the pyrochlore lattice connects exactly two β chains, but there are two inequivalent kinds of hexagon, as shown in <ref> (a). Neglecting constant terms, that is terms that do not distinguish different chain configurations {η_i}, the energy correction from a hexagon is proportional to the product of the Ising variables of the two chains i,j it connects
  
  
      δ E^(6)_ij   = const + J_chη_iη_j, 
      δ E^'(6)_ij   = const + J_ch'/2η_iη_j,
  
  
  where we use J_ch and J_ch' to distinguish the two ways in which two nearest-neighbor chains can be connected by a hexagon. The factor of 1/2 here is needed since there are 2L equivalent hexagons connecting to chains i and j in the primed direction (indicated in pink in <ref>). The interaction energy between two β chains i,j is then
  
      E^(int)_ij = L J_ch^(') η_i η_j,
  
  where 2L is the number of sites on the chain.
  
  The β chains embedded in the pyrochlore lattice thus form an anisotropic triangular lattice Ising model which can be described using an effective Hamiltonian
  
      ℋ_eff = J_ch∑_ijη_i η_j + J_ch' ∑_ij'η_i η_j,
  
  where the sum over ij runs over nearest neighbor pairs in the [112] and [1̅1̅2] directions and the sum over ij' indicates summation over nearest-neighbor pairs in the [110] direction. Note that here, we use “nearest-neighbor” to mean “connected by a hexagon” instead of referring to real-space distance since the β chains are farther apart in the [110] direction than in the other two directions. We indicate which chain coupling is generated by which linked cluster by matching colors in 
  <ref> (a) and (b). 
  
  The effective couplings J_ ch^(') are functions of the original pyrochlore exchange couplings J_λ, the field h and the mixing angle θ. The explicit functional form is different in each of the different chain phases. We thus denote the chain couplings in the Chain_λ phase by J_ ch^λ('):
  
      J_ ch^(') = 
              J_ ch^λ(')  in Chain_λ phase.
  
  We calculate the J_ ch^λ(') in sixth order RSPT (see app:rspt for details on the perturbative expansion). First consider the result in the Chain_Z phase
  
  
      J_ch^z =    - S J_z cos(ϕ)^2 J_x^2 J_y^2/J_z^4( J_x - J_y)^2/J_z^21/16Γ_z(h),
      
          J_ch^z' =    S J_z cos(ϕ)^4J_x^2 J_y^3/J_z^5( J_x/J_z + tan(ϕ)^2 ) 2/16Γ_z'(h),
  
  
  
  where S is the spin length and the J_λ are the exchange couplings in <ref>. The angle ϕ parameterizes the classical configuration of spins on the α chains such that θ-ϕ is the angle between the local field and the magnetic moments. ϕ thus depends on the external field h, θ as well as on the couplings J_x,J_z and ϕ→θ as h →∞. Γ_z and Γ_z' are dimensionless, positive factors controlling the asymptotic field dependence of the couplings. While full expressions are given in app:rspt, we note that the field dependence of Γ_z and Γ_z' is qualitatively different:
  
  
      Γ_z(h)∼ h^-2 as  h→∞,
      Γ_z'(h)∼ h^-4 as  h→∞,
  
  
  
  where the asymptotic exponents corresponds to the number of polarized sites on the hexagon generating the respective coupling [<ref> (a)]. This enables tuning the order-by-disorder strength and even the ground state by means of the magnetic field (see <ref> for details). 
  
  The respective expressions for J_ ch^x and J_ ch^x' are obtained from <ref> by interchanging J_x ↔ J_y and cosϕ↔sinϕ. It is also given explicitly in app:rspt.
  
  In the Chain_Y phase, the couplings in sixth order RSPT are given by
  
  
      J_ ch^y    = -1/64 S J_y sin(2ϕ)^2 J_x^2 J_z^2 (J_x - J_z)^2/J_y^6Γ_y(h), 
      
          J_ ch^y'   = 0,
  
  
  
  where Γ_y again is a dimensionless, positive, field dependent factor that vanishes as h^-2 as h →∞. For S=1/2, a finite contribution to J_ ch^y' even at higher orders is only generated when including a transverse term corresponding to cubic magnon terms in LSWT.
  
  The phase diagram of the effective model [<ref>] in terms of the effective couplings J_ch and J_ch' <cit.> is shown in <ref> (c). 
  There are three possible phases: if J_ ch is negative and -J_ ch>J_ ch' the system orders ferromagnetically. If the exchange parameter J_ ch is positive and larger than J_ ch', the system enters a `staggered' phase where there is ferromagnetic order along the primed direction (indicated in <ref> (b) in pink) with antiferromagnetic order in the perpendicular direction. Finally, if J_ ch' > J_ ch, the system is frustrated since when the system is ordered antiferromagnetically in the primed direction (indicated in <ref> (b) in pink), it is not possible to satisfy the non-primed bonds (indicated in <ref> (b) in green). Thus, in this `zigzag' phase, the ground state is highly degenerate, with each antiferromagnetically ordered row retaining an effective Ising degree of freedom. 
  Following the ferromagnetic order in the non-primed direction yields the eponymous zigzag shape, indicated in red in the sketch in <ref> (c).
  
  Since it is clear from <ref> and <ref> that J_ch is negative for all parameters, the chains will never order in a staggered configuration. This already yields a ferromagnetic ground state across the whole Chain_Y phase since there J_ch'=0 for S=1/2 at sixth order.
  In contrast, in the Chain_X and Chain_Z phases J_ch' is finite and can take either sign. The system is thus either in the ferromagnetic or zigzag phase, depending on sign and strength of J_ch' relative to J_ch.
  Focusing on the Chain_Z phase for simplicity, in terms of the original pyrochlore exchange couplings (J_λ), if J_y and J_x + J_z tan(ϕ)^2 have different signs, J_ch' is negative and we expect the chains to order ferromagnetically.
  In contrast, for J_x = J_y, J_ch vanishes while J_ch' is positive and finite, at least for J_x >0 and J_x < -J_z tan^2ϕ so we expect the chains to order in a zigzag configuration. 
  As indicated in <ref> (c), the phase boundary between the ferromagnetic and zigzag phases is obtained by equating J_ch' = -J_ch. This phase boundary as a function of the J_λ is also indicated in <ref> and <ref> as a solid line.
  
  Note that, while the leading order contribution to the effective chain couplings vanishes in certain regions of the phase diagram, we expect that higher-order terms will contribute also in those regions and eventually lift the classical degeneracy. However, since the degeneracy lifting in these cases is of higher order, it is strongly suppressed.
  
  
  
  
  
  
  
   §.§ Degeneracy lifting in the zigzag phase
  
  
  Dipolar-octupolar pyrochlores in a magnetic field undergo a remarkable dimensional evolution when considering different energy scales. On the largest energy scale, corresponding to the bare parameters of the Hamiltonian, the system appears fully three dimensional. In an intermediate regime (between J_ ch, J_ ch' and  J_z) the system undergoes a dimensional reduction such that it can be described as an ensemble of noninteracting one-dimensional chains. 
  At the lowest energy scales, one might expect the order by disorder mechanism to restore the three dimensional nature of the system.
  
  The zigzag phase, however, in some region of the phase diagram adds an additional stage to this dimensional evolution. Viewed in the pyrochlore lattice, left-over degrees of freedom are effectively two-dimensional. Magnetic moments on the β chains are ordered in stripe patters within the [hl0] plane, however the effective interactions between the planes are frustrated. Thus, the system has a ground state with a subextensive degeneracy scaling with the linear system size L, each plane carrying an effective Ising degree of freedom.
  
  This degeneracy is broken only at 10th order in perturbation theory, by a linked cluster indicated by the blue dashed line in <ref> (a). This cluster generates a next-nearest-neighbor coupling J_ ch” in the [001] direction as indicated in <ref> (b) and, equivalently, a ring exchange term. Note that the other possible next-nearest-neighbor coupling, while also generated at 10th order, is still frustrated in the zigzag phase and does not lift the degeneracy. As shown in <ref> (d), depending on the sign of the next-nearest-neighbor coupling J_ch”, there are two different zigzag states selected, either a `kinked' state, where ordered moments form a zigzag pattern in the [001] direction, or a `striped' state, where magnetic moments order ferromagnetically along one of the non-primed nearest-neighbor directions.
  
  The `ferromagnetic', `kinked' and `striped' configurations all have only an intensive number of degenerate states related by symmetries of the Hamiltonian [<ref>]. Hence, the system in the zigzag phase regains its full three dimensional correlations on an energy scale ∼ J_ ch” while on scales between J_ ch” and J_ ch, J_ ch' it features effectively two-dimensional degrees of freedom.
  
  We summarize the resulting hierarchy of energy scales in ground state selection in <ref>, where we also give rough numerical ranges in which the scales Δ and B lie for most values of exchange couplings if h ∼ J_z [cf. <ref> and <ref>].
  
  
  
  
  
  § ORDER BY DISORDER PHASE DIAGRAM
  
  
  
  
  
  
  The classical phase diagram of dipolar-octupolar pyrochlores in a moderate external field in the [110] direction, as discussed in <ref>, contains extended regions of chain phases in which the field generates effective one-dimensional degrees of freedom and which feature a subextensive degeneracy. Since the degenerate states are not related by symmetries, one expects the degeneracy to be lifted by quantum fluctuations. 
  In most of these chain phases (the Chain_λ phases), this degeneracy lifting, as established in <ref>, can be cast in the form of an anisotropic triangular-lattice Ising model with three effective chain coupling parameters. Motivated from the effective model, which was derived in leading-order real space perturbation theory (RSPT), we expect that on the largest energy scale, two different types of chain configuration will be chosen as a ground state. The chains will order either ferromagnetically, or in a zigzag state [<ref> (c)]. In the latter case, there is a further dimensional evolution to effectively two-dimensional degrees of freedom, which are planes perpendicular to the [001] direction, each carrying an effective Ising degree of freedom, noninteracting in leading-order RSPT. An effective interaction J_ ch”, breaking the residual degeneracy, is only generated at 10th order
  of RSPT. Depending on the sign of J_ ch”, one of = two zigzag configurations will be chosen, either a `kinked' or a `striped' state both of which then feature an intensive number of degenerate states related by symmetries [<ref> (d)].
  
  
  
  
  
   §.§ Phase diagram
  
  
  In this section, we turn to calculate the full classical degeneracy lifting in the Chain_λ phases, as a function of the exchange parameters J_λ. To corroborate the RSPT calculation, we calculate the lifting numerically using linear spin wave theory (LSWT). The calculation is performed using a semi-infinite slice of a pyrochlore lattice, that is a 8× 8 triangular lattice of infinite β chains with periodic boundary conditions in the directions orthogonal to the chains.
  
  The central quantity to compute ground state selection by quantum fluctuations is the zero point energy
  
      ℰ_0 ( {S⃗_i^(0}) = 1/π∑_ν∫_0^π/2 k ω_ν((k, -k, 0)),
  
  where ω_ν(q⃗) are the magnon frequencies of the system around the ground state {S⃗_i^(0)}. Details of the numerical computation can be found in app:lswt.
  We compute the zero point energy for all possible ground state configurations of the effective nearest-neighbor model, which are the ferrromagnetic state, the staggered state and 12 different zigzag states (those are all zigzag states on a 8× 8 triangular lattice not related by lattice symmetries). 
  Assuming that the effective model is a valid low-energy description of the system, this will yield the correct ground state phase diagram and serve as a quantitative test of the effective couplings calculated in RSPT. 
  To test the validity of the effective model, we also compute the zero point energies for 1000 randomly selected chain configurations and compare them to the predictions of the effective model, with three couplings fitted using LSWT (see <ref> for details). We find that the effective model indeed gives an accurate estimate of the energies, and therefore that searching for ground states only among the possible ground states of the effective model is justified. 
  
  The resulting (quantum) order-by-disorder phase diagram for the Chain_λ (λ = x,y,z) phases is shown in <ref> and <ref> for θ=0.25π and θ=0.33π respectively. The selected configurations are those expected from the effective Ising model. Remarkably, the phase boundaries obtained from LSWT also match very well with those obtained directly from RSPT at least for moderate transverse couplings. 
  Going beyond the RSPT calculation, LSWT also reveals the splitting of the zigzag state into the striped and kinked configurations.
  There are also parts of the phase diagram where the classical degeneracy is not lifted within LSWT. Most notably, this occurs along the two lines J_x = 0 and J_y = 0 for which we prove explicitly that ω(q⃗) are identical for all classical ground states within LSWT in app:lswt.
  
  To shed light on the relevant energy scales in the ground state selection by quantum fluctuations, we also show in <ref> and <ref> the gap between the ferromagnetic and the zigzag configurations Δ as well as the bandwidth of the zigzag states B (that is also the gap between the striped and the kinked state).
  It is evident that the energy scales on which the classical degeneracy is lifted are very small and have a clear hierarchy as illustrated in <ref>. 
  In particular, for the compound Nd_2Zr_2O_7, using the most recent estimate of the exchange couplings<cit.> and a field of about 1 Tesla yields Δ≈ 10^-8 meV. In this example, quantum OBD has thus no experimental relevance in ground state selection. 
  
  To investigate the influence of finite temperature on OBD, we employ two approaches (see <ref> for details). First, we compute the leading order contribution of the classical low-temperature expansion of the free energy and show that it is identical for all classical ground states in the Chain_λ phases. Second, we compute the free energy of the noninteracting magnon gas at finite temperature and show that below a crossover temperature T_ co set by the spin wave gap, the energy scales Δ and B are barely modified from their zero temperature values. For Nd_2Zr_2O_7, we estimate T_ co≈ 0.23 K (see <ref>).
  
  Our findings are consistent with the fact that disordered chain configurations were observed in neutron scattering experiments on this compound <cit.> at low temperatures.
  
  
  
   §.§ Validity of effective model
  
  
  In order to further assess the validity of the effective model, we compare it to the full spectrum of zero point energies.
  To this end, we compute the effective chain couplings J_ ch and J_ ch' explicitly using the zero point energies obtained from LSWT using
  
      J_ ch   = 1/4( ℰ_0^ FM - ℰ_0^ Staggered), 
      
          J_ ch'    = 1/2( ℰ_0^ FM + ℰ_0^ Staggered), 
      
          J_ ch”   = 1/2( ℰ_0^ Kinked - ℰ_0^ Striped).
  
  We then compare the zero point energies ℰ_0({S⃗_i^(0)}) of all 1014 classical ground state configurations that we consider (ferromagnet, staggered, 12 zigzag states and 1000 randomly chosen configurations) to the estimate from the effective nearest-neighbor model using the effective couplings obtained from Eq. <ref>. 
  
  Within the Chain_Z and Chain_X phases, the estimates agree remarkably well with the energies computed from LSWT, with residues between the two on the order of the next nearest-neighbor coupling J_ ch” [see <ref> in app:lswt].
  Restricting ourselves to the 12 zigzag states, their splitting is modeled by the fitted J_ ch” up to 10^-12 J_z.
  
  In the Chain_Y phase, we expect the effective model to be still valid as derived in RSPT. At high fields, both LSWT and RSPT predict dominance of the ferromagnetic configuration for all exchange couplings in that phase. However, at low fields LSWT for some part of the phase diagram predicts a zigzag ground state, in contrast with RSPT. 
  The disagreement can be rationalized by considering that LSWT is unable to capture the hard-core constraint of magnons for S=1/2 and this constraints leads the coupling J_ ch' to vanish at sixth order in RSPT.
  Meanwhile, in RSPT we neglect higher orders which could lead to a finite J_ ch' even for S=1/2.
  
  
  
  
  
  
  
  § FIELD TUNABILITY
  
  
  
  
  As is qualitatively clear from <ref>, the effective nearest-neighbor couplings depend on the external magnetic field h. 
  Since the effective interaction between chains is generated by fluctuations on the α chains, the external field can directly control the strength of these fluctuations and hence be used to manipulate the OBD effect. This is clearly seen in <ref>, where the difference of the zero point energies of the ferromagnetic state and the zigzag states is shown as a function of field for two sets of exchange parameters.
  
  Furthermore, the ground state that is selected by quantum fluctuations can itself be tuned by means of the field. This is the case because selection of the ferromagnetic state or the zigzag states is determined by a competition of the effective couplings in the two inequivalent nearest-neighbor directions of the lattice. Inspecting the two leading-order contributions shown as green and pink hexagons in <ref> (a), we see that the two hexagons have two and four sites located on α chains respectively. Hence, as is reflected in the different behaviour of Γ_z and Γ_z' in <ref>, the couplings scale differently with the field and hence the competition between J_ ch and J_ ch' changes as a function of field. As an example of this, we show in the top panel of <ref> how the energy difference between the ferromagnetic state and the zizag states (computed using LSWT) changes sign as a function of field.
  
  
  
  
  
  
  § EXPERIMENTAL RELEVANCE
  
  
  
  In this Section we discuss the implications of our results for the dipolar-octupolar pyrochlore compounds Nd_2Zr_2O_7 and Ce_2Sn_2O_7.
  
  
  
   §.§ Disordered dipolar chain state in Nd2Zr2O7
  
  
  First, a transition from an AIAO ordered phase to a disordered chain phase induced by a field has been observed in Nd_2Zr_2O_7 <cit.>. Such a transition is expected for -J_x < J_z < -J_x/3 or -J_z < J_x < -J_z/3, which is consistent with the exchange parameters of this compound as estimated from neutron scattering, for which the most recent estimate (taken from Ref. Xu2019) is 
  
      J_x ≈ 0.1 meV,   J_z ≈ -J_x/2,   J_y ≈ 0.15 J_x.
  
  
  For these exchange parameters and assuming a field of 1 Tesla, LSWT predicts a ferromagnetic order of the chains that is however separated from the zigzag state by a gap of only Δ≈ 10^-8 meV. At temperatures above T_ co≈ 0.23 K, this is even further suppressed (see discussion in <ref>). Hence, even at the lowest experimentally accessible temperatures one would not expect to see ordering of the chains driven by quantum fluctuations, which is consistent with the experimental observation of a disordered chain state.
  
  This makes Nd_2Zr_2O_7 in a moderate [110] field an excellent platform for the study of one-dimensional quantum XYZ chains.
  
  
  
   §.§ Octupolar quantum spin liquid in Ce2Sn2O7
  
  
  
  
  Second, recent experimental evidence for an octupolar U(1) quantum spin liquid (QSL) has been reported in the dipolar-octupolar pyrochlore Ce_2Sn_2O_7 <cit.>. This implies that in this compound the octupolar exchange is positive and large <cit.>
  
      J_y > max(J_x, J_z, -3J_x, -3J_z).
  
  
  Classically, this would result in an ice-like phase of octupolar moments SI_Y at zero-field. This phase has the same extensive degeneracy as spin-ice. At finite fields, the system would be in the Chain_Y^* phase in which the extensive degeneracy is broken down to a subextensive degeneracy with one Ising degree of freedom carried by each α and each β chain [<ref>]. At a finite critical field, the system would then undergo a second-order transition to the Chain_Y phase, which also features a subextensive degeneracy but with only half as many Ising degrees of freedom, located on the β chains. 
  This is shown in <ref>, where we plot the (classical) octupolar staggered moment [<ref>] as a function of the external field for exchange parameters as used to model Ce_2Sn_2O_7 in Ref. Sibille2020.
  The second order transition at finite field is absent if the ground state is a dipolar spin ice state SI_X and SI_Z since those can be directly connected to the respective chain phases Chain_X and Chain_Z.
  
  In a real compound, where S=1/2, quantum fluctuations might be expected to be strong, the story laid out above gets slightly modified. The zero field octupolar ice phase SI_Y will be replaced by a U(1) QSL that prevails for some finite range of external field<cit.>. At some critical field h_c^(1) this is followed by a confinement transition to a quantum version of the Chain_Y^* phase. At a larger, second critical field h_c^(2), there will then be a transition to the `simple' chain phase Chain_Y. 
  
  The transition from the QSL to Chain_Y^* can be captured  by means of degenerate perturbation theory around the low field limit. As we show in app:quantum-chain-y-star, for S=1/2 there is a diagonal perturbation at fourth order in the external field that generates an effective antiferromagnetic (with respect to the local basis) interaction within the α chains. 
  This interaction favours a set of states with one-dimensional degrees of freedom on each α and β chain, corresponding precisely to the Chain_Y^* phase found in our classical analysis.
  This indicates that the Chain_Y^* phase exists in the quantum limit, and should appear in octupolar quantum spin ices subject to moderate [110] fields.
  
  The double transition from a U(1) QSL to the Chain_Y^* phase and finally to the Chain_Y phase as a function of a [110] field appears to be unique to the case of dominantly octupolar exchange and can hence serve as another experimental test for compounds suspected to realize such physics like Ce_2Sn_2O_7.
  Additionally, the existence of the Chain_Y^* phase is interesting in and of itself. It serves as an example of how uniform external fields acting on noncollinear magnetic states can lead to the development of complex degrees of freedom. While the number of degrees of freedom grows only with the square of the linear system size L^2, their interaction should not be expected to be modeled by a local two-dimensional effective model since α and β chains interpenetrate in a nontrivial way.
   
  
  
  
  
  
  § CONCLUSION
  
  
  In conclusion, we have mapped out the classical phase diagram of dipolar-octupolar pyrochlores in a [110] field in detail and studied quantum order-by-disorder (OBD) in the `chain' phases, in which the classical ground state degeneracy can be parametrized by a subextensive number of effective Ising degrees of freedom, each describing the magnetic moment of a one-dimensional chain. We have focused our study on a subset of these phases, where all degrees of freedom are carried by a set of parallel chains. We show that OBD in this case can be modelled remarkably well by a simple effective anisotropic triangular lattice Ising model that we derive in real space perturbation theory (RSPT) and corroborate using numerical linear spin wave theory (LSWT). 
  
  Ground state selection by quantum fluctuations in dipolar-octupolar pyrochlores is a multi-step process governed by a hierarchy of energy scales. First, competition of two inequivalent nearest-neighbor chain interactions drives selection of either a ferromagnetic ordering of the chains or a `zigzag' order. In the latter case, the chains are ordered antiferromagnetically in one direction while the nearest-neighbor coupling in the other two directions is frustrated, leaving a subextensive degeneracy scaling with the linear system size. This left-over degeneracy is finally broken by an effective next-nearest-neighbor coupling, albeit on a minuscule energy scale.
  Taken together, the system undergoes a highly nontrivial dimensional evolution when considering effective degrees of freedom on different energy scales.
  Comparing the results of the effective model with LSWT, the two methods mostly agree remarkably well even quantitatively, with phase boundaries between ferromagnetic and zigzag phases almost indiscernible in the dipolar chain phases Chain_X, Chain_Z for moderate transverse couplings. 
  
  
  Highlighting connection to different experiments, we have also discussed that the energy scales on which quantum fluctuations break the classical degeneracies are minuscule for parameter values relevant to Nd_2Zr_2O_7. This is consistent with the experimental observation of a disordered chain state in this compound <cit.> and makes this material a promising platform to study pure quantum XYZ chains. 
  Finally, we study the field-breakdown of the octupolar quantum spin liquid (QSL) state recently observed in Ce_2Sn_2O_7. We predict that for a field in the [110] direction, the transition to the large-field chain state Chain_Y happens in two steps, via an intermediate phase that we call Chain_Y^*. This sequence of transitions is unique for the case of dominantly octupolar exchange and can hence be used as another test for the nature of Ce_2Sn_2O_7.
  
  
  Acknowledgements:
  This work was in part supported by the Deutsche Forschungsgemeinschaft  under grants SFB 1143 (project-id 247310070) and the cluster of excellence ct.qmat (EXC 2147, project-id 390858490). 
  
  
  
  
  
  
  
  
  § DERIVATION OF LEADING-ORDER CORRECTION IN REAL SPACE PERTURBATION THEORY
  
  
  
  
   §.§ Setup
  
  
  In this appendix, we derive the expression for the leading-order effective chain couplings presented in <ref> using real space perturbation theory. We start from a Hamiltonian of the form
  
      ℋ = ∑_ijS⃗_i ·J⃗·S⃗_j - ∑_i h⃗_i ·S⃗_i,
  
  from which <ref> is obtained for the choice J⃗ = diag(J_x, J_y, J_z) and h⃗_i = μ_B g_z H⃗· z_i (sinθ,  0,  cosθ),
  with H∥ (1,1,0). 
  
  We transform to a local basis, u_i, v_i, w_i chosen such that S⃗_i^(0) = S w_i is a classical ground state.
  The Hamiltonian is then separated an unperturbed part, ℋ_0 and four transverse perturbations
  
      ℋ = ℋ_0 + ∑_ij V_ij^(1) + V_ij^(2) + V_ij^(3) +V_ij^(4),
  
  where
  
  
      ℋ_0    = E_0 + ∑_i B_i δS_i, 
  
      E_0    = S^2/2 ∑_i ( ∑_j∈nn(i) w_i ·J⃗ ·w_j ) - S w_i ·h⃗_i, 
  
      B⃗_i    = -S ∑_j∈nn(i) J⃗ ·w_j + h⃗_i,
  
  
  and
  
  
      V_ij^(1)   = c⃗_i^* ·J⃗·c⃗_j^*  S^+_i S^+_j + h.c.,
          
      
          V_ij^(2)   = c⃗_i^* ·J⃗·c⃗_j  S^+_i S^-_j + h.c.,
          
      
          V_ij^(3)   =  w_i ·J⃗· w_j δ S_i δ S_j,
          
      
          V_i^(4)   = c⃗_i ·J⃗· w_i  S^-_i δ S_j +  w_i ·J⃗·c⃗_j δ S_i S^-_j + h.c.,
  
  
  
  where δ S_i = S - S_i^w and c⃗_i = 12( u_i + i  v_i).  
  
  
  Going forward, we will first calculate the effective couplings in the Chain_Z phase (see <ref>), with the corresponding results for the Chain_X phase being obtained by the same procedure. 
  Finally, the effective couplings in the Chain_Y phase will be computed in the last part of this appendix.
  
  For a ground state in the Chain_Z phase (see <ref>) we have
  
  
      w_i =
          (σ_i^αsinϕ, 0, σ_i^αcosϕ)    i∈α chain
      (0, 0, σ_i^βη_c )    i∈β chain,
  
  and choose
  
      c⃗_i = 1/2(σ_i^αcosϕ, i, -σ_i^αsinϕ)    i∈α chain
      (1, iσ_i^βη_c, 0 )    i∈β chain,
  
  
  
  where η_c is the Ising variable of the chain that the site i belongs to
  
  
      σ_i^α   =
          
              1    if  i ∈ℒ_0 
      
              -1    if  i ∈ℒ_3
          ,
      σ_i^β   =
          
              1    if  i ∈ℒ_1 
      
              -1    if  i ∈ℒ_2
          .
  
  
  
  The angle ϕ above parameterizes the classical configuration of spins on the α chains such that θ-ϕ is the angle between the local field and the magnetic moments. ϕ thus depends on the external field h, θ as well as on the couplings J_x,J_z and ϕ→θ as h →∞.
  
  The perturbations V_ij^(1…4) hence depend parametrically on the chain configuration {η_i} via the basis <ref>. Standard perturbation theory yields a configuration dependent energy correction Δ E({η_i }). As we will see below, the leading order contribution to Δ E depends only on products of nearest-neighbor chains such that
  
      Δ E = ∑_ijδ E_ij =  J_ ch∑_ijη_iη_j + J_ ch' ∑_ij'η_iη_j,
  
  which yields the effective Ising model discussed in <ref>.
  
  
  
  
   §.§ Perturbation theory
  
  
  Since the energy corrections must stay extensive, any term that contributes to perturbation theory must correspond to a linked cluster of the pyrochlore lattice, where each edge represents one or more applications of the perturbation V = V^(1) + V^(2) + V^(3) + V^(4). The energy correction corresponding to a linked cluster with p edges is <cit.>
  
      δ E^(p)  = ∑_{ψ_i}⟨0| V |ψ_1⟩⟨ψ_1| V |ψ_2⟩…⟨ψ_p-1| V |0⟩/(E_0 - E_ψ_1) … (E_0 - E_ψ_p-1),
  
  where E_0 is the classical ground state energy and E_ψ_i is the unperturbed energy of the intermediate state |ψ_i⟩. From the above we can read off some simple rules:
  
      
    * Any linked cluster corresponds to a sum of terms, which are all possible combinations of applications of the V_ij^(1…4) to the links such that we begin and end in the vacuum of spin flips |0⟩.
      
    * Any such sequence of perturbations has to begin and end with either V^(1) or V^(4).
      
    * Any such sequence must involve an even number of applications of V^(4).
  
  Since we are interested only in the leading order contribution to the energy correction we do not have to consider the perturbations V^(3) and V^(4) since the operator δ S_i does not modify the spin configuration. Thus for any sequence of applications of perturbations involving V^(3) or V^(4), there will be a nonvanishing lower order term that contributes.
  
  The matrix elements of the perturbations are evaluated straightforwardly using <ref> and <ref>. If i,j are both on a β chain
  
  
      V_ij^(1)   = 1/4( J_x + J_y) S_i^+ S_j^+ + h.c., 
      
          V_ij^(2)   = 1/4( J_x - J_y) S_i^+ S_j^- + h.c..
  
  In contrast, if i,j are both on an α chain 
  
      V_ij^(1)   = -1/4(cosϕ^2 J_x + J_y + sinϕ^2 J_z) S_i^+ S_j^+ + h.c., 
      
          V_ij^(2)   = 1/4(-cosϕ^2 J_x + J_y - sinϕ^2 J_z) S_i^+ S_j^- + h.c..
  
  Finally, if i is on a β chain and j is on an α chain
  
      V_ij^(1)   = 1/4(σ^α_j cosϕ J_x - σ^β_iη_c J_y) S_i^+ S_j^+ + h.c., 
      
          V_ij^(2)   = 1/4(σ^α_j cosϕ J_x + σ^β_iη_c J_y) S_i^+ S_j^- + h.c..
  
  
  
  
  Similarly, the factors in the denominator of <ref> are given by
  
      E_0 - E_ψ   = ⟨ψ|E_0 - H_0|ψ⟩
         = -∑_i B_i ⟨ψ|δ S_i |ψ⟩,
  
  where if i is on a β chain:
  
  
      B_i = 2 S J_z,
  
  and if i is on an α chain:
  
      B_i  = 2 S (sin(ϕ)^2 J_x + cos(ϕ)^2 J_z) + cos(ϕ-θ)  h.
  
  
  
  
  
  
   §.§ Degeneracy lifting
  
  
  We are interested in the lowest order quantum correction that lifts the classical degeneracy. For that, the energy correction must be proportional to the product of as least two chains η_i η_j, which is only possible if the corresponding cluster connects two chains. 
  
  
  
  The two smallest clusters that connect two chains are shown in <ref> (a). They both have two links and their smallest contribution is of 4th order. Inspecting <ref>, we see that any term that is proportional to η_c also has a factor σ^β_i. Hence, the two clusters marked in red and blue will cancel out exactly. This will be true for any cluster that has a `mirrored' cluster that connects to a site on the same chain but on a different sublattice.
  Hence, only clusters of nontrivial topology, i.e. those that have noncontractible loops, will contribute to the effective chain couplings. The smallest such cluster on the pyrochlore lattice is a hexagon and there are two inequivalent ways a hexagon can connect two chains, shown in <ref> (a). Note that for the same reasons, there will be no nonvanishing contribution from perturbation theory, that is proportional to only a single chain variable η_i. Any such term must be proportional to a respective σ^β_i and has a mirrored term that gives the same contribution but with the opposite sign.
  
  Evaluating the energy contributions of the two types of hexagons, we have to consider all possible distributions of V_ij^(1) and V_ij^(2). However, not all of those contribute to the chain coupling. Consider for example the contribution from the cluster depicted in <ref> (b), which is [using <ref>]
  
      δ E^(6)∝   (cosϕ J_x + η_1 J_y ) (cosϕ J_x - η_1 J_y ) 
         (cosϕ J_x + η_2 J_y ) (cosϕ J_x - η_2 J_y )(J_x - J_y)^2
  
  and thus the dependence on η_i drops out completely.
  
  Keeping only those configurations which do depend on the chain configuration {η_i}, there are four inequivalent possibilities, shown in <ref> (c). Each of them have different possible sequences of intermediate states {ψ_i}, which must be enumerated to evaluate the denominator of <ref>. We show one possible such sequence, for the first configuration in <ref> (c), in <ref> (d), where the excited states are shown as dots on the sites.
  
  Finally, the effective coupling J_ ch is given as the sum of the four inequivalent contributions depicted in <ref> (c)
  
      J_ch =   4(2S)^6/4^6cos(ϕ)^2 J_x^2 J_y^2 ( (J_x - J_y)^2 (Γ_1 - Γ_2)  
         + (J_x + J_y)^2 (Γ_3 - Γ_4) ),
  
  where we have defined factors Γ_i, stemming from the denominator of <ref>
  
      Γ_config := -∑_{ψ_i}1/(E_0 - E_ψ_1) … (E_0 - E_ψ_p-1)
  
  Note that since all factors in the denominator are negative and there are an odd number of factors for each closed loop on the pyrochlore lattice, Γ_ config>0.
  
  Evaluating the Γ factors, one finds that Γ_3 = Γ_4 and hence 
  
  
      J_ch^z =    - S J_z cos(ϕ)^2 J_x^2 J_y^2/J_z^4( J_x - J_y)^2/J_z^21/16Γ_z(h),
  
  where Γ_z(h) := S^5 J_z^5 (Γ_2 - Γ_1) is defined as a dimensionless factor controlling the asymptotic field dependence of the couplings and we have introduced a superscript `z' to denote that the above result is true in the Chain_Z phase. 
  
  
  
  Computing J_ch' follows exactly the same procedure. Since the sites on the hexagon, that is the leading-order contribution to J_ch' are on different sublattices, also configurations of perturbation on the hexagon that contribute will be different. The four inequivalent configurations that contribute to J_ch^z', together with their contributions are shown in <ref>. Summing those up yields 
  
      J_ch^z' =    S J_z cosϕ^4J_x^2 J_y^3/J_z^5( J_x/J_z + tanϕ^2 ) 2/16Γ_z'(h).
  
  
  
  This completes the derivation of the results stated in the main text. For completeness, we give the Γ-factors explicitly:
  
  
      Γ_z(h)/(SJ_z)^5 =   2/B^β_z (B^α_z + B^β_z)^4 
              + 2/(B^β_z)^2 (B^α_z + B^β_z)^3
         + 1/(B^β_z)^3 (B^α_z + B^β_z)^2, 
      Γ'_z(h)/(SJ_z)^5 =   4/B^β_z (B^α_z + B^β_z)^4,
  
  
  
  where B_z^α and B_z^β are the local fields on the α and β chains respectively, as defined in <ref>.
  Note that since the local fields B^α_λ and B^β_λ are always positive, so are Γ_z and Γ_z'
  
  In the Chain_X phase, the ground state as defined in <ref> will be slightly different 
  
  
      w_i =
          (σ_i^αsinϕ, 0, σ_i^αcosϕ)    i∈α chain
      (σ_i^βη_c, 0, 0 )    i∈β chain,
  
  and
  
      c⃗_i = 1/2(σ_i^αcosϕ, i, -σ_i^αsinϕ)    i∈α chain
      (0, iσ_i^βη_c, -1 )    i∈β chain.
  
  
  
  Note the sign change in the definition of c⃗_i on β chains that is necessary so the basis stays right-handed. Since J⃗ is diagonal, the above change can be accounted for by simply changing J_x ↔ J_z and cosϕ↔sinϕ in <ref> and hence
  
  
      J_ch^x =    - S J_x sin(ϕ)^2 J_z^2 J_y^2/J_x^4( J_z - J_y)^2/J_x^21/16Γ_x(h), 
      
          J_ch^x' =    S J_x sin(ϕ)^4J_x^2 J_y^3/J_x^5( J_z/J_x + cotan ϕ^2 ) 2/16Γ_x'(h).
  
  
  
  where Γ_x^(') is obtained from Γ_z^(') [<ref>] by substituting B_z^α and B_z^β by B_x^α and B_x^β respectively.
  
  
  
  
   §.§ Chain-Y phase
  
  
  We now turn to calculate the effective couplings in the Chain_Y phase, which is done mostly by the same procedure as used for the Chain_X/Z phases. 
  As before, we start by defining a basis where the w_i are given by the classical ground state configuration
  
  
      w_i =
          (σ_i^αsinϕ, 0, σ_i^αcosϕ)    i∈α chain
      (0, σ_i^βη_c, 0 )    i∈β chain,
  
  and the c⃗_i complete the basis
  
      c⃗_i = 1/2(σ_i^αcosϕ, i, -σ_i^αsinϕ)    i∈α chain
      (i, 0, σ_i^βη_c )    i∈β chain.
  
  
  
  Substituting the above into <ref> yields if i,j are both on a β chain
  
  
      V_ij^(1)   = -1/4( J_x + J_z) S_i^+ S_j^+ + h.c.,
      
          V_ij^(2)   = 1/4( J_x - J_z) S_i^+ S_j^- + h.c.,
  
  in contrast, if i,j are both on an α-chain
  
      V_ij^(1)   = -1/4(cosϕ^2 J_x + J_y + sinϕ^2 J_z) S_i^+ S_j^+ + h.c., 
      
          V_ij^(2)   = -1/4(cosϕ^2 J_x - J_y + sinϕ^2 J_z) S_i^+ S_j^- + h.c.,
  
  and finally if i is on a β-chain and j is on an α-chain
  
      V_ij^(1)   = -1/4σ^α_j ( icosϕ  J_x + σ^β_iη_c sinϕ  J_z) S_i^+ S_j^+ + h.c., 
      
          V_ij^(2)   = -1/4σ^α_j ( icosϕ  J_x + σ^β_iη_c sinϕ  J_z) S_i^+ S_j^- + h.c..
  
  
  
  Note that in the last case, the matrix elements are now complex.
  
  
  
    §.§.§ Effective coupling Jch-y
  
  
  We now turn to identify the degeneracy-breaking perturbations. In the case of the Chain_Z/X phases those are given by linked clusters of nontrivial topology. Note that the argument for this relied solely on the fact that any Ising chain variable η_c comes also with a factor σ^β_i, which is clearly still the case. Hence we can again restrict the discussion to the two kinds of hexagons as long as we are interested only in leading-order corrections. 
  
  
  
  As before, not all distributions of the perturbations V_ij^(1) and V_ij^(2) yield an energy correction that is proportional to the two Ising variables η_1, η_2 on the hexagon. For that, first consider the case of J_ ch^y, which is generated by the hexagon marked in green in <ref> (a). In this case, proportionality to a chain variable η_1 comes from a product of two perturbations connecting that chain to a polarized site, where the factor σ^β_i is different for both. Inspecting <ref>, we see that the two perturbations have to act differently on the two β chain sites in order for the resulting energy correction to be proportional to η_1. Note that this rule already excludes the possibility of having V^(1)_ij act on the β chain bonds [see <ref> (a)] and hence J_ ch∝ (J_x-J_z)^2 in accordance with the LSWT result (see <ref> and <ref>).
  Fixing the perturbations acting on β chains to spin-flip hopping, there are four diagrams contributing to J_ ch^y shown in <ref> (b). Summing up their contributions yields
  
      J_ ch^y    = -1/64 S J_y sin(2ϕ)^2 J_x^2 J_z^2 (J_x - J_z)^2/J_y^6Γ_y(h), 
      Γ_y    = S^5 J_y^5 ( 2Γ_2^Y + 2Γ_3^Y - 2Γ_1^Y - Γ_4^Y),
  
  where the Γ^Y_i are computed as before by enumerating all possible intermediate states and Γ_y is given by
  
      Γ_y/(SJ_y)^5 =   2/(B_y^β)^3 (B_y^α + B_y^β)^2 - 2/(B_y^α + 2B_y^β)^3(B_y^α + B_y^β)^2
         + 4/(B_y^β)^2 (B_y^α + B_y^β)^3 - 4/(B_y^α + 2B_y^β)^2 (B_y^α + B_y^β)^3
         + 2/B_y^β (B_y^α + B_y^β)^4 - 2/(B_y^α + 2B_y^β)(B_y^α + B_y^β)^4,
  
  where B_y^α and B_y^β are the local fields on the α  and β chains respectively in the Chain_Y phase, obtained from <ref> by replacing J_z by J_y. Note that since the local fields are always positive, it is clear that Γ_y>0 and hence the effective chain coupling J_ ch^y<0 is ferromagnetic for all fields and exchange parameters.
  
  
  
    §.§.§ Effective coupling Jch-y'
  
  
  We now turn to the second effective exchange coupling, J_ ch^y', which is generated by the hexagon marked in pink in <ref> (a). On the level of RSPT diagrams, this amounts to a different distribution of polarized and unpolarized sites on the hexagon as indicated in <ref> (c). Turning to identify the degeneracy breaking perturbations, we straightforwardly arrive at a similar rule as in the case of J_ ch^y. However, since the two perturbations connecting the site on the β chain now connect to the same site, and because σ_j^α appears just as an overall sign in the matrix elements [<ref>], the two perturbations now have to act in the same way (i.e. either both with S^+ of both with S^-) on the site on the β chains as indicated in <ref> (c). 
  However, this is not possible at sixth order since there can be no double occupation of excitations on any site. 
  For the same reason, for S = 1/2 any diagram using only V_ij^(1), V_ij^(2) cannot yield a finite contribution to J_ ch^y', even at higher orders.
  Since in the Chain_Y phase V_ij^(3)=0 if i on a β chain and j on an α chain (or the other way around), a finite contribution to J_ ch' can only arise from a perturbation including V_ij^(4). Such a term occurs at seventh order in RSPT and competes with contributions from higher orders to J_ ch, even when the sixth order contribution to J_ ch vanishes.
  
  Since in dipolar-octupolar pyrochlores the pseudospin always forms a doublet, we conclude that J_ ch^y' is strongly suppressed for these compounds. Taken together with the fact that J_ ch^y < 0, we conclude that quantum fluctuations should mostly or even always select a ferromagnetic chain configuration in the Chain_Y phase. 
  LSWT (see <ref> and <ref>) reproduces the overall dominance of ferromagnetic ordering at high fields, but predicts a zigzag ground state for some exchange parameters at low fields. We ascribe this to the fact that LSWT does not capture the hardcore nature of magnons at S=1/2 since V_ij^(4) corresponds to cubic magnon terms in spin wave theory.
  
  Notwithstanding the ambiguity of the two methods at low fields, at higher fields both LSWT and RSPT calculations agree that ferromagnetic chain configurations will dominate the Chain_Y phase [see <ref> (b) and <ref> (b)].
  
  
  
  
  
  
  
  § LINEAR SPIN WAVE THEORY
  
  
  In this appendix, we provide details about the linear spin wave theory calculations presented in the main text.
  We start from a system consisting of L unit cells U_l (l=0… L-1), repeated along one dimension. Each unit cell has M sites at positions l a⃗ + r⃗_i (i=0… M-1). 
  The general nearest-neighbor exchange Hamiltonian of such a system reads
  
      ℋ =   1/2∑_l ∑_i,j∈ U_lS⃗_i 𝒥_ijS⃗_j
              + ∑_l ∑_i∈ U_l∑_j∈ U_l+1S⃗_i 𝒦_ijS⃗_j 
         - ∑_l ∑_i ∈ U_lh⃗_i ·S⃗_i,
  
  where 𝒥_ij and  𝒦_ij parametrize the exchange within and between unit cells respectively.
  
  To derive the linear spin wave Hamiltonian, we first write the spins in a local basis u_i, v_i, w_i, chosen such that in the classical ground state all spins align with w_i
  
      S⃗_i = S_i^uu⃗_i + S_i^vv⃗_i + S_i^ww⃗_i,
  
  and define raising and lowering operators via the (linear) Holstein-Primakoff transformation
  
  
      S_i^u   ≈√(S/2)(a_i + a_i^†), 
      
          S_i^v   ≈√(S/2)(a_i - a_i^†), 
      
          S_i^w   = S - a_i^† a_i.
  
  
  Then, we perform a one-dimensional Fourier transform on sublattice i
  
      a^†_i(q⃗) = 1/√(L)∑_l exp(-i q⃗· (l a⃗  + r⃗_i)),
  
  where the momentum takes values q⃗ = 2πa⃗ k/ a⃗^2 L (k=0… L-1). Using the above, we arrive at
  
      ℋ =    E_0 + ℋ_1 + ℋ_2,
  
  where E_0 is the classical ground state energy, ℋ_1 contains terms linear in the operators a_i, a_i^† and vanishes if S⃗_i^(0) = S w_i is a classical ground state. Finally, ℋ_2 is given by
  
  
      ℋ_2    = 1/2∑_q⃗A⃗(q⃗)^†𝒳(q⃗) A⃗(q⃗), 
      A⃗^†(q⃗)     = (a_0^†(q⃗), …, a_M-1^†(q⃗), a_0(-q⃗), …), 
      𝒳(q⃗)    = S (X^(11)(q⃗)     X^(12)(q⃗) 
       X^(21)(q⃗)     X^(22)(q⃗)),
  
  
  
  
      X_ij^(11)(q⃗)    = e^-i q⃗·r⃗_ij(
                  c⃗_i 𝒥_ijc⃗_j^*
                  + e^-i q⃗·a⃗c⃗_i 𝒦_ijc⃗_j^*
                  + e^i q⃗·a⃗c⃗_j^* 𝒦_jic⃗_i 
         =e^-i q⃗·r⃗_ij(
                      - δ_ij(
                          ∑_l=0… M-1(
                              w⃗_i 𝒥_ilw⃗_l
                              + w⃗_i 𝒦_ilw⃗_l
                              + w⃗_l 𝒦_liw⃗_i
                          )
                          - 1/Sh⃗_i ·w⃗_i
                      )
               ), 
      
          X_ij^(12)(q⃗)    = e^-i q⃗·r⃗_ij(
              c⃗_i 𝒥_ijc⃗_j
              + e^-i q⃗·a⃗c⃗_i 𝒦_ijc⃗_j
              + e^i q⃗·a⃗c⃗_j 𝒦_jic⃗_i
              ), 
      
          X_ij^(21)(q⃗)    = e^-i q⃗·r⃗_ij(
              c⃗_i^* 𝒥_ijc⃗_j^*
              + e^-i q⃗·a⃗c⃗_i^* 𝒦_ijc⃗_j^*
              + e^i q⃗·a⃗c⃗_j^* 𝒦_jic⃗_i^*  ), 
      
          X_ij^(22)(q⃗)    = e^-i q⃗·r⃗_ij(
                  c⃗_i^* 𝒥_ijc⃗_j
                  + e^-i q⃗·a⃗c⃗_i^* 𝒦_ijc⃗_j
                  + e^i q⃗·a⃗c⃗_j 𝒦_jic⃗_i^* 
         =e^-i q⃗·r⃗_ij(
                      - δ_ij(
                          ∑_l=0… M-1(
                              w⃗_i 𝒥_ilw⃗_l
                              + w⃗_i 𝒦_ilw⃗_l
                              + w⃗_l 𝒦_liw⃗_i
                          )
                          - 1/Sh⃗_i ·w⃗_i
                      )
               ),
  
  
  
  
  where the X^(nm)(q⃗) are M× M matrices. The sum goes over q⃗ = 2πa⃗ k/ a⃗^2 L (k=0… L-1) and a⃗ is the translation vector between neighboring unit cells and r⃗_ij = r⃗_j - r⃗_i is the difference vector between two sites within the unit same cell.
  
  <ref> is quadratic and thus can be diagonalized using a Bogoliubov transformation. Explicitly, the magnon dispersion is obtained by diagonalizing σ𝒳(q⃗) where
  
      σ = (𝕀_M× M    0 
       0     -𝕀_M× M).
  
  The matrix σ𝒳(q⃗) has eigenvalues ω_ν(q⃗), -ω_ν^*(q⃗), ν=0… M-1.
  
  The magnon spectrum ω_ν(q⃗) depends on the classical ground state {S⃗_i^(0} used in the expansion above. The ground state selection by quantum fluctuations can then be quantified using the zero point energy
  
      ℰ_0 ( {S⃗_i^(0})    = lim_L→∞1/L⟨0| H_2 |0⟩
         = a⃗/4π∑_ν∫_0^2π/a⃗ k ω_ν(k a⃗ / a⃗).
  
  
  
  
   §.§ Details of the numerical calculations
  
  
  Here, we provide some details of the computation of the numerical phase diagram in <ref> and <ref>. For a discussion of the different phases, we refer to the main text [<ref>].
  
  As mentioned already in the main text, the LSWT theory calculations are performed on a semi-infinite “tube” of the pyrochlore lattice, that is an effective W× W triangular lattice of β chains. The infinite direction is thus [11̅0] [a⃗ = (2, -2, 0)] and the cluster has M = 4 W^2 sites per unit cell. 
  The zero point energy is then given by
  
      ℰ_0 ( {S⃗_i^(0})    = 1/π∑_ν∫_0^π/2 k ω_ν((k, -k, 0)),
  
  where we used that the dispersion is symmetric around k=0.
  
  For a fixed set of exchange couplings J⃗ and field h⃗, we diagonalize the resulting 2M× 2M matrix numerically for 50 values of q⃗ and compute the integral of the dispersion numerically using the standard Simpson's quadrature rule. Since relative differences between different zero point energies are very small, we verified that using 50 points, the integral is converged up to double precision. For W=8, computing the zero point energy of a single classical ground state for fixed parameters takes about 10 seconds.
  
  The number of distinct configurations of the Ising variables {η_i} grows exponentially with the system size W. The last system size for which it is feasible to enumerate all possible configurations is W=4, for which there are 2^W^2=65 536. This can be reduced to 674 inequivalent chain configurations which are not related by symmetries. Computing the zero point energies [<ref>] for all these states in linear spin wave theory yields the phase boundary between the ferromagnetic and the zigzag state, as well as the gap Δ between the two. 
  However, at this system size, LSWT is not able to correctly capture the splitting between zigzag states. This is because the unit cell of the two chosen ground states, the kinked and the striped states is of size W=4. If we would simulate the effective triangular lattice Ising model directly, a system of the size of that unit cell would be sufficient to capture the selection process.
  
  However, since we actually consider the whole pyrochlore lattice and wish to capture the
  fluctuations which generate the interactions between the chains, we need a cluster which is at least twice as large to suppress the influence of virtual processes which cross the entire cluster via the periodic boundary conditions.
  
  
  In order to obtain the fully resolved phase diagram shown in <ref> and <ref>, we use a system with size W=8. For each set of parameters, we compute the zero point energy of all possible ground states of the effective model, that are the ferromagnet, the staggered state, and all inequivalent zigzag states. The latter is possible because the total number of zigzag states grows exponentially only in the linear system size, hence there are 2^W=256 zigzag states in total. The number can be further reduced to 12 zigzag states not related by symmetry. To summarize, there are 14 possible inequivalent ground states of the effective model for W=8.
  
  
  
   
  
  To benchmark the validity of the effective triangular lattice Ising model we calculate the zero-point energy of 1000 random states. We show the results of this benchmarking in  <ref>. In (a) we compare the zero point energies obtained by LSWT to those computed from an effective model with effective couplings fitted using three states [<ref>]. In (b), we compare the effective chain couplings obtained from 6th order RSPT (J_ ch^(')(6)) to those computed from LSWT (J_ ch^(') SW).
  
  Finally, to assess finite size effects we also simulate a system with W=16 and verify that the numerical values of Δ and B are independent of the size.
  
  
  
  
   §.§ Survival of full classical degeneracy when Jx=0 or Jy=0 in the dipolar chain phases
  
  
  In this appendix, we prove the survival of the full classical degeneracy, within LSWT, in the singular cases of J_x=0 or J_y=0. To this end, we relate the flipping of single β chain in the system to a unitary transformation on the matrix σ𝒳(q⃗). This already implies that the spectrum and hence the zero point energy of the system does not depend on the chain configuration {η_i} and hence the full classical degeneracy is preserved even in the presence of quantum zero-point fluctuations.
  
  In the dipolar chain phase i.e. in either Chain_X or Chain_Z, all classical ground states are periodic in the [11̅0] direction with a unit cell of four sites per β chain in the system. Two sites of those are on the β chain and the other two are on α chains.
  We now consider the case of a general (nonzero) matrix element of 𝒳(q⃗) as given in <ref>. Flipping a single chain variable η_c → -η_c amounts to a change in a maximum of 18 matrix elements in each of the four submatrices X^(mn) of 𝒳(q⃗). These matrix elements can take two forms.
  
  In the following, we assume for simplicity to be in the Chain_Z phase, but the same argument applies in the Chain_X phase. First, if i,j are on the same β chain, the matrix elements are given by
  
  
      X^(11)_ij(q⃗)    = X^(22)_ij(q⃗) = cos(π k /L) (J_x - J_y), 
      
          X^(12)_ij(q⃗)    = X^(21)_ij(q⃗) = cos(π k /L) (J_x + J_y),
  
  
  where we used the definition of c⃗_i from <ref>. The above is independent of the Ising variable η_c of the chain and hence there are only 16 matrix elements left that might change as the chain is flipped. 
  The other sixteen matrix elements are all of the form 
  
  
      X^(11)_ij(q⃗)    = X^(22)_ij(q⃗) 
         = 1/2exp(± iq⃗·a⃗/4)
              ( σ_j^α J_x cosϕ + σ_i^βη_c J_y ), 
      
          X^(12)_ij(q⃗)    = X^(21)_ij(q⃗) 
         = -1/2exp(± iq⃗·a⃗/4)
              ( σ_j^α J_x cosϕ - σ_i^βη_c J_y ),
  
  
  
  where i is a site on a β chain, j is a site on an α chain and σ_i^α, σ_j^β are defined in <ref>. The sign of the site-dependent phase ± i q⃗·a⃗/4 depends on the specific sites i and j, but this explicit dependence does not matter for the argument presented here.
  For J_y = 0, it is clear that the spin wave matrix 𝒳(q⃗) does not depend on the Ising variable η_c and hence there is no ground state selection. For J_x=0 the argument is less obvious but the same is still true. In that case, all matrix elements that change when η_c → -η_c are proportional to η_c. Thus, flipping the chain is equivalent to a simple unitary transformation on σ𝒳(q⃗) with the transformation matrix given by 
  
      U = diag(1, …, -1, -1, …, -1,-1, … 1),
  
  where the -1 entries appear at the positions of the two sites on the flipped β chain.
  
  
  
  
  
  § PERTURBATIVE LIMIT OF THE CHAIN-Y* PHASE 
  
  
  In this appendix, we show that the Chain_Y^* phase can be derived also in degenerate perturbation theory, from the limit of the octupolar ice phase SI_Y. 
  
  
  For simplicity, we set the transverse coupling to zero and study the Hamiltonian
  
  
      ℋ   = ℋ_0 + V, 
      ℋ_0    = J_y ∑_ij S_i^y S_j^y, 
       
          V    = - h ∑_i∈ℒ_0(S_i^+ + S_i^- ) + h ∑_i∈ℒ_3(S_i^+ + S_i^- ).
  
  
  
  The ground states of ℋ_0 are given by the octupolar spin-ice configurations, where ∑_i∈ tetS_i^y = 0 on every tetrahedron. 
  
  The leading order correction appears at second order and consists of the perturbation V_i acting on the same site twice
  
      Δ E_i^(2)   = ⟨ψ_0| V_i ( E - H_0 )^-1 V_i |ψ_0⟩
         = 2Sh^2/-2S J_y = -h^2/J_y.
  
  The above is just a constant correction to the unperturbed ground state energy E_0, that is it does not discriminate between ground states ψ_0, leaving the degeneracy unchanged. 
  
  This changes at fourth order, where the only nonvanishing contribution consists of application of the perturbation V_i to two neighboring sites. Note that the perturbation acts only on the α chains. In the following, we assume that S=1/2, however qualitatively the result is the same for other values of S. The energy correction is then given by
  
      Δ E_i^(4)   = 4 ×(2S)^2 h^4/-(2SJ_y)^2 S J_y (3 + σ_i σ_j)
         = const + h^4/2J_y^3σ_i σ_j,
  
  where σ_i=±1 denotes the configuration of the octupolar moment on site i in the ground state ψ_0. The above expression now does discriminate between ground states ψ_0, favoring those states with a staggered configuration of octupolar moments on the α chains.
  
  The ice rule ∑_i ∈ tet S^y_i=0 then forces all β chains into a staggered configuration 
  This leaves an Ising degree of freedom on each of the α chains and each of the β chains, reproducing exactly the degeneracy structure of the classical Chain_Y^* phase.
  
  
  
  
  The above result can be easily understood. First, since any valid sequence flips two spins twice, the numerator is given by (2S)^2h^4. 
  A possible sequence of intermediate states is shown in <ref> with contributions from each tetrahedron to the energy relative to the ground state indicated in red. Since the ground state is not allowed as an intermediate state, there are exactly three other possible sequences which all have the same sequence of intermediate energies. 
  The first and last state have one of the two moments S_i^y or S_j^y flipped, and hence have an energy of E-E_0=2SJ_y.
  The second intermediate state has both moments flipped and thus violates the spin-ice rule on three tetrahedra. Two of those include just one of the two sites and hence contribute also a term 2SJ_y. The tetrahedron that shares the two sites however contributes an energy depending on the ground state configuration. If the two flipped sites where initially oriented antiparallel, then the spin ice rule is still fulfilled on this tetrahedron if both are flipped, while in the case of the two sites being parallel, the spin ice rule is now maximally violated. Altogether the energy of the second intermediate state can thus be written as E-E_0=SJ_y(3+σ_iσ_j).
  
  
  
  Adding transverse exchange (J_x, J_z) to <ref> will allow the appearance of a ring exchange term in perturbation theory which will favour a U(1) QSL <cit.>.
  The competition between this ring exchange and the diagonal term Δ E_ij^(4) will determine the low field phase boundary between the QSL and the Chain_Y^∗ phase.
  Since the transition from QSL to Chain_Y^∗ does not require the condensation of spinons, but instead confines their dynamics it is a type of confinement transition.
  This confinement transition should precede the Anderson-Higgs type transition from spinon condensation predicted in gauge Mean Field Theory calculations <cit.>, which generally requires h∼ J_y in order to close the spinon gap.
  
  
  
  
  § EFFECT OF FINITE TEMPERATURE
  
  
  In this appendix, we study the effect of finite temperature of ground state selection. First we consider the low-temperature expansion of the classical free energy around different ground states and show that the first order correction is independent of the ground state. Second, we consider the free energy of the magnons at finite temperature.
  
  We show that this leads to the same ground state selection
  for all temperatures.
  For temperatures below a crossover temperature T_ co, set by the spin wave gap, the strength of ground state selection is close to the zero temperature value. For  T>T_ co the difference in free energy between different chain configurations is heavily suppressed. For parameters of interest (close to those estimated for Nd_2Zr_2O_7) this crossover temperature is given by T_ co≈ 0.2 J_z (see <ref>).
  
  
  
   §.§ Low temperature expansion
  
  
  For the low temperature expansion, we choose as a starting point the same general one-dimensional Hamiltonian as for linear spin wave theory (app:lswt).
  
  That is a system consisting of L unit cells U_l (l=0… L-1), repeated along one dimension. Each unit cell has M sites at positions l a⃗ + r⃗_i (i=0… M-1) and the Hamiltonian reads
  
      ℋ =   1/2∑_l ∑_i,j∈ U_lS⃗_i 𝒥_ijS⃗_j
              + ∑_l ∑_i∈ U_l∑_j∈ U_l+1S⃗_i 𝒦_ijS⃗_j 
         - ∑_l ∑_i ∈ U_lh⃗_i ·S⃗_i,
  
  We again proceed by writing the spin variables ins a local basis u⃗_i, v⃗_i, w⃗_i chosen such that in the (classical) ground state all spins align with w⃗_i. We can hence parametrize fluctuations around the ground state as
  
      S⃗_i    = √(S) δ u_i u⃗_i  +
              √(S) δ v_i v⃗_i +
              √(S^2 - S (δ u_i)^2 - S (δ v_i)^2 )w⃗_i 
         ≈√(S) δ u_i u⃗_i +  
              √(S) δ v_i v⃗_i +
              (S - 12(δ u_i)^2 - 12(δ v_i)^2) w⃗_i
  
  substituting the above into <ref> up to quadratic order in the fluctuations δ u_i, δ v_i yields
  
  
      ℋ = E_0 + ℋ_1 + ℋ_2
  
  where E_0 is the ground state energy, ℋ_1 is linear in the fluctuations δ u_i, δ v_i and vanishes if S⃗^(0)_i = S w⃗_i is a classical ground state, and ℋ_2 is given by
  
      ℋ_2    = 1/2∑_q⃗ũ⃗̃(-q⃗)^T 𝔛(q⃗)  ũ⃗̃(q⃗), 
      ũ⃗̃(q⃗)^T    = (δ u_1(q), δ u_2(q), …, δ u_M(q), δ v_1(q), …).
  
  
  
  Here, the matrix 𝔛(q⃗) is related to the spin wave Hamiltonian 𝒳(q⃗) <ref> derived in app:lswt by a unitary transformation
  
  
      𝔛(q⃗)    = F^†𝒳(q⃗) F 
       
          F    = 1/√(2)(𝕀   -i𝕀
      𝕀    -i𝕀)
  
  
  
  
  The classical partition function can be calculated explicitly
  
      Z    =  e^-βℋ = (2π)^-ML∫(∏_ν=1^2Mũ_ν) e^-β E_0 
              e^-βℋ_2
         = e^-β E_0∏_q⃗(β𝔛(q⃗) )^12
  
  and the free energy up to linear order in temperature by
  
      F    = -T log Z 
         = E_0 + T/2∑_q⃗log𝔛(q⃗) - TML log T.
  
  
  Now, note that flipping a chain in all cases amounts to flipping w⃗_i on two sites. To keep the basis right handed, one then also has to flip either u⃗_i or v⃗_i on those sites. 
  In the expression for the partition function (<ref>) however, assuming we flipped  u⃗_i, one can revert this change by the transformation of variables δ u_i → -δ u_i. Since the sign of an even number of variables is flipped, this leaves the phase space volume element (∏_ν=1^2Mũ_ν) invariant. Hence, the partition function Z as well as the free energy F are invariant under flipping any single chain variable η_c up to second order in deviations from the ground state configuration. Generally the above argument implies that for a set of (locally) collinear ground states as considered here in all of the chain phases, there is no order-by-disorder when taking the purely classical limit.
  
  
  
   §.§ Free energy of magnons at finite temperature
  
  
  
  
  The second approach treating the finite temperature is considering the gas of non-interacting magnons. For that, we start from the spin wave Hamiltonian in <ref>, which after a Bogoliubov transformation takes the form
  
      ℋ_2 = ∑_q⃗, νω_ν(q⃗) ( n_ν(q⃗) + 1/2)
  
  where n_ν(q⃗) = η_ν^†(q⃗) η_ν(q⃗) is the magnon number operator and the magnon dispersion ω_ν(q⃗) is given by the eigenvalues of σ𝒳(q⃗). Note that these eigenvalue come in pairs ±ω_ν(q⃗) (see app:lswt for details). 
  
  The free energy of this system at temperature T is then
  
      F_2    = -T log( e^-βℋ_2) 
         = ∑_q⃗, ν Tlog(1-e^-βω_ν(q⃗))+1/2ω_ν(q⃗)
  
  where the second term in the sum yields the zero point energy and the first term is a finite-temperature correction. For zero temperature, the above is equal to the zero point energy while for large temperature it becomes independent of the classical ground state configuration {S⃗^(0)_j}. Particularly, for all temperatures the same ground state is favored by fluctuation of the free Bose gas but the gap in free energy vanishes as T→∞. As an example, we show the free energy gap Δ between the ferromagnetic chain configuration and the zigzag band as a function of temperature, for exchange parameters as estimated for Nd_2Zr_2O_7, in <ref>.
  
  The vanishing order-by-disorder in the high temperature limit can be understood by expanding the exponential inside the logarithm in <ref>
  
      F_2 =
         ∑_q⃗, ν Tlog(βω_ν(q⃗) - 1/2(βω_ν(q⃗))^2 + (βω_ν(q⃗))^3)
          
          
          
          +1/2ω_ν(q⃗) 
         = ∑_q⃗, ν Tlog[ω_ν(q⃗)] - Tlog(T) + βω_ν(q⃗)^2.
  
  Note that the Eigenvalues of σ𝒳(q⃗) are exactly the pairs ±ω_ν(q⃗) and the number of eigenvalues M is even such that σ𝒳(q⃗) is positive. Hence, we can write the leading order term of F_2 as
  
      T ∑_q⃗, νlog[ω_ν(q⃗)]    = T/2∑_q⃗logσ𝒳(q⃗).
  
  Now, using that σ = 1 and that the quantum spin wave matrix 𝒳(q⃗) is related to its classical counterpart 𝔛(q⃗) (defined in <ref>) by a unitary transformation [<ref>], we arrive at
  
      T ∑_q⃗, νlog[ω_ν(q⃗)] = T/2∑_q⃗log𝔛(q⃗).
  
  At infinite temperature, we hence recover the purely classical result where, as discussed in the previous section of this appendix, there is no order by disorder.
  
  We note that the discussion of finite-temperature OBD as presented here is not comprehensive as it does not incorporate the possible presence of domain walls at finite temperature. However, since in the presence of anisotropies domain wall excitations, just as spin wave excitations, are gapped, we expect our results to be robust in the low-temperature limit.
  
  
  
  
  
  