Theory of magnetostriction for multipolar quantum spin ice in pyrochlore materials

Multipolar magnetism is an emerging ﬁeld of quantum materials research. The building blocks of multipolar phenomena are magnetic ions with a non-Kramers doublet, where the orbital and spin degrees of freedom are inextricably intertwined, leading to unusual spin-orbital entangled states. The detection of such subtle forms of matter has, however, been difﬁcult due to a limited number of appropriate experimental tools. In this work, motivated by a recent magnetostriction experiment on Pr 2 Zr 2 O 7 , we theoretically investigate how multipolar quantum spin ice, an elusive three-dimensional quantum spin liquid, can be detected using magnetostriction, by examining the characteristic signatures of its magnetic-ﬁeld descendent multipolar kagome ice phase, as well as that of the neighboring multipolar ordered phases in the pyrochlore materials. We provide theoretical results based on classical and / or quantum studies of non-Kramers and Kramers magnetic ions, and contrast the behaviors of distinct phases in both systems. Our work paves an important avenue for future identiﬁcation of exotic ground states in multipolar systems.


I. INTRODUCTION
The development of a robust understanding of emergent phenomena in strongly correlated quantum systems [1,2] ultimately requires deep insight into the many-body ground state, and the excitations it can support.This historically successful paradigm has been challenged in recent times by two prominent examples: quantum spin liquids (QSLs) and multipolar ordered states (MPOs).QSLs, which are long-range entangled correlated paramagnets, support deconfined fractionalized excitations (spinons) that couple to an emergent gauge field [3,4].QSLs arise from a variety of mechanisms, from geometrical frustration [5,6] to anisotropic bond-dependent interactions [7][8][9].The lack of magnetic ordering presents an obvious challenge as to its detection with conventional probes, and despite efforts from neutron scattering [10,11], a true smoking-gun signature has proven to be elusive.Analogously, MPOs also defy detection by conventional probes, despite falling under the purview of spontaneous symmetry breaking.These ordered states, which arise from spin-orbit coupling and crystalline electric fields (CEFs) placing restrictions on localized electron orbitals' shapes, do not possess just a simple dipolar moment.Instead, they support nontrivial charge and magnetic density distributions (described by higher rank multipolar moments [12,13]), which fail to directly couple to neutrons and other probes of ordering, and have been appropriately named "hidden orders." The central question that remains for both of these phenomena is the following: How can the existence and properties of QSLs' and MPOs' interacting many-body ground states be examined if they shy away from conventional probing tools?Some answers to this question have recently been found, where novel elastic-based techniques in multipolar heavy fermion systems seem to indicate the onset of MPOs [14][15][16][17].Motivated by such experiments, we ask whether the experimentally reticent QSLs can be exposed if they arise from interacting multipolar moments.
In this work, we propose that magnetostriction (length change under an external magnetic field) provides a sharp and distinct signature of quantum spin ice formed from non-Kramers multipolar moments.The current work is directly motivated by a recent experiment on the quantum spin ice candidate material, Pr 2 Zr 2 O 7 [35].In order to validate this proposal, we contrast its difference by presenting the distinctive length change behaviors of possible ordered states in Kramers and non-Kramers ions.In doing so, we establish a comprehensive theory of magnetostriction for a number of possible emergent phases in both non-Kramers and Kramers pyrochlore systems.The findings are based on corroborating classical analysis and exact diagonalization of quantum models on the pyrochlore lattice.Our theoretical results provide a means to identify the existence of a QSL's many-body ground state, by focusing on the signatures of its field-induced descendent multipolar kagome ice phase, and to distinguish the various ordered phases in the pyrochlore family.Our work lays the foundations upon which targeted experimental investigations can be conducted on QSLs, and provides a new direction of inquiry in the field of multipolar magnetism.

II. PSEUDOSPIN-1/2 MODEL OF PYROCHLORE MATERIALS
In the R 2 M 2 O 7 oxide family, the local moments arise from the f electrons of the R 3+ rare-earth-metal ions.Importantly, the surrounding cage of O 2− ions subject the f electrons to a local D 3d crystalline electric field (CEF), which splits the J multiplet of isolated R 3+ ions to yield low-lying Kramers or non-Kramers doublet ground states, depending on the nature of the R 3+ ions.In Yb 2 Ti 2 O 7 [24,[36][37][38][39], for example, the J = 7/2 multiplet is split to yield Kramers ground states that support conventional magnetic dipole moments, which can be efficiently represented by the pseudospin-1/2 operator S = J.In the intriguing candidate quantum spin-ice material, Pr 2 Zr 2 O 7 , the J = 4 degenerate manifold is partially lifted to yield non-Kramers (doublet) ground states of an even number of f electrons.As a consequence, these support, in addition to a conventional magnetic dipole moment, more exotic timereversal even quadrupolar moments.The multipolar moments can be efficiently represented by the pseudospin components, S x = J x J z , S y = J y J z , and S z = J z , where the overline indicates a symmetrized product [26].In the local frame of each sublattice (as described in Appendix A), members of the R 2 M 2 O 7 family obey the following generic nearest-neighbor pseudospin-1/2 model [26,27,40], where i j is a sum over nearest neighbor sites i and j on the pyrochlore lattice, γ i j and ζ i j = −γ * i j are unimodular complex numbers listed in Appendix A. J zz > 0 is an antiferromagnetic interaction that gives rise to the celebrated two-in, two-out spin ice rule.Because of the quadrupolar nature of S x,y in the non-Kramers variety, J z± = 0 by time-reversal symmetry.We note that, in anticipation of the discussion to follow, the exchange coupling constants are conventionally taken to be independent of applied magnetic fields, as the groundstate doublet in Pr 2 Zr 2 O 7 is well separated from its excited states [10].
Figure 1 presents the T = 0 classical phase diagram associated with Eq. (1) with J z± = 0.This classical non-Kramers phase diagram provides a zoo of possible phases: a classi-FIG.1. Classical phase diagram of Eq. ( 1) for non-Kramers ions (J z± = 0).The depicted phases are multipolar spin ice (MSI) of T 1g , coplanar antiferroquadrupolar (cAFQ L ) of T 2g , a second coplanar antiferroquadrupolar (cAFQ L ) of T 1g , and ferroquadrupolar (FQ L ) of E g symmetry.Here we use the subscript L to indicate orderings in the local basis (Appendixes A and N).cAFQ L and cAFQ L are related to each other by a local C 4z rotation on each sublattice.cal two-in, two-out multipolar spin ice (MSI) phase of T 1g symmetry, a coplanar antiferroquadrupolar (cAFQ L ) phase with T 2g symmetry, another coplanar antiferro-quadrupolar (cAFQ L ) of T 1g symmetry, and a ferroquadrupolar (FQ L ) of E g symmetry [41].Here we use the subscript L to indicate orderings in the local basis.
Although we have so far discussed classical multipolar phases, it is highly suggestive from parton mean-field theory (known as gauge mean field theory, gMFT [30,31]) studies that quantum phases are the descendants of these parent classical phases.Indeed, in gMFT, the classical SI phase gets promoted to a U(1) quantum spin-liquid phase, which is characterized by the existence of deconfined bosonic spinons (magnetic monopoles, in the spin-ice literature [42]) coupled to a U(1) gauge field [28].The other classically ordered phases get promoted to Higgs phases in gMFT, where the bosonic monopole condenses thus eliminating the emergent gauge field.In the current work, instead of using gMFT, we examine the quantum model using exact diagonalization, which indeed confirms the relevant phase diagram (Appendix H).

III. ELASTIC STRAIN COUPLING TO LOCAL MOMENTS
In this section, we examine the coupling of the local R 3+ moments to the elastic normal modes.The cubic nature of the underlying Bravais lattice constrains (by O h point group) the elastic energy to be of the form, where the crystal's deformation is described by the components of the strain tensor ik , and c i j is the elastic modulus tensor describing the stiffness of the crystal.Here c B is the bulk modulus, B ≡ xx + yy + zz is the volume expansion of the crystal, ν ≡ (2 zz − xx − yy )/ √ 3 and μ ≡ ( xx − yy ) are cubic normal mode lattice strains.Because of the sublattice nature of pyrochlore lattice, we specify that the elastic strain tensors, magnetic fields, as well as local moments, written in the local basis of a given sublattice α possess a sublattice index, i.e., S x,y,z α , h x,y,z α , α i j .In the global basis {(1, 0, 0), (0, 1, 0), (0, 0, 1)}, we write the same quantities without the sublattice index.In coupling the elastic strain to the localized moments, we enforce the local D 3d point group symmetry of the surrounding CEF locally.We present in Appendixes B and C the relationships between local and global quantities and their transformations under D 3d .
First, since the non-Kramers XY pseudospin components are time-reversal even electric quadrupolar moments, they can couple linearly to elastic strain as where we have introduced Einstein summation notation for α.
We explicitly denote the non-Kramers case by the subscript NK, and k 1,2 are phenomenological coupling constants.The Z pseudospin component contains the time-reversal odd magnetic dipole moment and can only couple to the elastic strain in the presence of a time-reversal breaking external magnetic field, h, to yield xx + α yy , (4) where g 1,2,3,4 are coupling constants, and we again employ Einstein summation notation for α.Equation ( 4) is common to both non-Kramers and Kramers ions.

IV. MAGNETOSTRICTION BEHAVIORS OF NON-KRAMERS PYROCHLORE MATERIALS
Under an applied magnetic field, the magnetic dipole moment couples at linear order and to the quadrupolar moments at quadratic order, where t sums over all up tetrahedra and S μ t,(α) is the S μ moment on sublattice α of tetrahedron t.The second term in Eq. ( 5) is a quadratic-in-h coupling to the quadrupolar moments, which is perturbatively weak as compared to the magnetic field coupling to the dipole moments; we include small phenomenological quantities δ 12, to represent this diminutive nature.The strength of the quadratic background in the length change depends on δ 1,2 .Indeed, the value of δ 1,2 depends inversely on the gap ( ) between the ground state and excited states, δ 1,2 ∼ 1/ .In the context of Pr 2 Zr 2 O 7 , the gap is relatively large ( ≈ 9.5 meV [10]), thus physically justifying the minuscule magnitude of δ 1,2 .In other pyrochlore materials, such as Tb 2 Ti 2 O 7 , this gap is almost an order of magnitude smaller ( ≈ 1.4 meV [44]), and consequently the quadratic background to the magnetostriction is more dominant [45,46].Applying the strategy described inAppendix D, we derive the following length change expressions along the = (1, 1, 1) and where the subscript and superscript in ( L L ) refer to the magnetic field and length change directions, respectively.We note that the pseudospin operators are taken to be understood as their expectation value with respect to the ground state.
We note that we have redefined the couplings in Eqs. ( 6) and (7) from Eqs. (3) and (4) for brevity, i.e., g 1 ≡ g 1 FIG. 2. Length change, L L , under applied [111] magnetic field, h for spin-ice phase J ± = 0.02J zz , J ±± = 0.05J zz (a) along the (1,1,1) direction and (b) along the (1,1,0) direction.For an infinitesimal field, the classical spin ice enters into the degenerate Kagome spin-ice phase (denoted by yellow shaded region).The quantum spin ice is, however, stable for a small window of magnetic field strengths [43]; we schematically denote this region by the orange shaded region (the size of the region is amplified for ease of viewing).For large fields, the system resides in a fully polarized state (indigo-blue shaded region).The green, blue, and red curves [squares] denote the length change arising from the XY pseudospin (quadrupolar), Z pseudospin (dipole), and combined contributions, respectively, from classical [32 site exact diagonalization] studies.The dipole plot line is (slightly) purposely shifted from the total length change plot line to more easily visualize the individual contributions.The classical (1,1,0) length change possesses two behaviors in the kagome ice phase, due to the degeneracy of the kagome ice manifold being reflected in the length change [Eq.( 7)].The average over the degenerate branches (red dashed line) matches up with the exact diagonalization (ED) result for the (1,1,0) direction.The values of the chosen lattice-pseudospin couplings are presented in Appendix F. ordering of the XY local moments results in vanishing length contributions from the quadrupolar moments.This scenario occurs for when J ±± = 0 and only J zz , J ± > 0. To have nonvanishing contributions from the quadrupolar moments, the = (1, 1, 1) length change clearly requires the assistance of J ±± = 0 to give a nonuniformity (or even canting) to the ordering on each sublattice.For = (1, 1, 0), a ferrolike ordering still yields a finite length change contribution from the quadrupolar moments.

A. Unique magnetostriction signature of multipolar spin ice
We present in Fig. 2 the unique magnetostriction behavior for the MSI phase along the (1,1,1) and (1,1,0) directions under an applied field along the [111] direction.The solid lines are obtained from a classical computation, while the squares are obtained from 32-site exact diagonalization (ED) of the quantum model.We describe in Appendixes G and H the procedure for these respective techniques.The depicted shaded regions can be understood as a battle between two energy scales: spin exchanges and magnetic field.For small fields, the classical system enters into a kagome ice (KI) phase, where the spin on sublattice 0 is fully polarized, while the remaining three spins conspire to satisfy the overall twoin, two-out ice rules of the exchange terms over an entire tetrahedron [47][48][49].For the quantum model, the U(1) QSL survives for small window of magnetic field strength [43] (depicted by an amplified orange shaded region in Fig. 2, for ease of viewing), until it enters into the quantum kagome ice phase (described below).In the classical KI phase, the extensive degeneracy of the classical MSI phase partially remains within each kagome layer [50,51].This can be easily noticed on a given tetrahedron where three possible states satisfy the twoin, two-out ice rules: ↓, ↓, ↑}.In the quantum limit, the ground state is that of a superposition over this degenerate manifold and is named quantum kagome ice, QKI (when discussed in the context of the simplified X X Z model, where J ±± = 0) [52,53].Finally, in the large-field limit, all the pseudospins are polarized: This state corresponds to the pseudospins on sublattice-0 (sublattice-1,2,3) pointing out of (into) a given tetrahedron.Interestingly, an island of XY quadrupolar ordering develops for intermediate fields during the transition between the KI and fully polarized state, where the XY orderings are not identical on the sublattices.This results in a sharp discontinuity in the XY length change, accompanied by a "flip" in the z component on sublattice 1 to S z < 0 (or any of the sublattices that have S z = 1/2), and an alignment into the fully polarized state: We present in Appendix I the behavior of the local pseudospin configurations on each sublattice under the [111] magnetic field.
Both classical total length changes (the sum of the quadrupolar and dipolar contributions) in Fig. 2 possess a sharp discontinuity which originates from the transition from KI to the fully polarized phase.Furthermore, both direc-tions possess an additional quadratic-in-h scaling behavior (imposed on top of a linear-in-h scaling from the dipole moments).This arises from the ∼h 2 coupling of the magnetic field to the quadrupolar moments in Eq. (5).A crucial difference between the two length change directions is mainly with the dipole contributions to the respective length change direction.First, in the fully polarized limit, (1,1,1) and (1,1,0) directions have opposite signs in their total length change.Second, and more interestingly, the (1,1,0) direction has two possible (dipole) length behaviors from the classical computation in the KI phase, while there is a unique behavior for the (1,1,1) direction.This is a result of the degeneracy of the KI phase.For the (1,1,1) direction, all three degenerate KI states (KI a , KI b , KI c ) give the same expressions for the length change when these configurations are inserted into Eq.( 6).For the (1,1,0) direction, however, the length-change expression of Eq. ( 7) is sensitive to which of the three KI degenerate states is chosen.In particular, (KI a , KI b ) have the same length-change behavior when their respective configurations are inserted in to Eq. ( 7), but (KI c ) has a different length-change expression.This difference can be traced to the dipole terms unique to the (1,1,0) direction in Eq. (7).Depending on which KI degenerate solution is chosen, we can thus get one of the two branches as depicted in Fig. 2(b), and subsequently the sign of the dipole contribution for the (1,1,0) direction can flip (or be retained) in the fully polarized limit.In a realistic system, it is possible that one may obtain an average over the three possible KI configurations.Interestingly, this "averaged" behavior is precisely what the ED computation of the quantum model finds, corroborating the idea that the quantum ground state can be thought as the superposition of three degenerate configurations.All of these behaviors thus provide a sharp signature for the existence of a MSI phase at h = 0.
The ED results match well with the classical solution's (1,1,1) direction, especially in the region of KI and the fully polarized state.Although the quadrupolar peak is broadened out as compared to its classical counterpart, we attribute this to possible finite-size effects of ED, quantum fluctuations, and the challenge of extracting the symmetry-broken order parameter.We provide a detailed explanation of the latter point in Appendix H.For the (1,1,0) direction, more care needs to be taken to understand its apparent difference with its classical counterpart.In particular, the obtained ED ground state is nondegenerate, with the z-spin expectation value on each of sublattice 1, 2, 3 being − 1 6 .This appears to suggest that the quantum ground state is an equal superposition over all the three degenerate classical kagome ice states (on a given tetrahedron), resulting in the single possibility for the dipole length change.For completeness, and to enable comparison with experiments, we present the length change under a [110] magnetic field in Appendix J. Comparing the relatively smooth magnetostriction features for the [110] field to the peak structure for [111] field highlights the strong anisotropy and selection rules of magnetostriction.

B. Length change behaviors of non-Kramers multipolar ordered phases
In order to clearly distinguish SI signatures from the other magnetically (or quadrupolar) ordered ground states, it is pertinent to consider the magnetostriction behavior of the symmetry-broken ordered phases.It is beneficial to first write, in terms of the classical multipolar order parameters (using Appendix N), both the interacting pseudospin model, and the = (1, 1, 1) magnetostriction expression (Eq.6),

L L
[111] (1,1,1),NK where we drop the gerade subscript for the order parameters (described in Fig. 1) for brevity, and collect the constants under , and . Figure 3 depicts the magnetostriction behavior for = (1, 1, 1) of the various multipolar ordered phases discussed earlier; the vertical dashed lines indicate jump discontinuous behaviors in the ordering.Specifically, both the cAFQ L and cAFQ L have jump discontinuous behaviors which correspond to S z (0) becoming fully polarized.The cAFQ L also has the distinction of having a finite length change in the absence of an external field.This is apparent from Eq. ( 9), where the cAFQ L order parameter is present even for h = 0.The FQ L± states do not possess any nonanalytic behavior in their length changes.Indeed, the local moments undergo a smooth and gradual change into the fully polarized state.We note that the FQ L+ and FQ L− behaviors are related by a local C 4z rotation of the pseudospins, where the ± denote J ±± > 0 and J ±± < 0, respectively; these different parameter choices of the same phase highlight the independence of the qualitative features of the magnetostriction on the precise value of the exchange couplings.
Because of these mentioned characteristics, each of the MPOs have their own distinct signature that allows each of them to be identified individually, as well as be distinguished from MSI.The FQ L± states are the easiest to identify, as they possess a smooth change in the length change; this gradual change is not present in any of MSI nor the other MPOs.cAFQ L can also be distinguished as it holds the honor of being the only non-Kramers phase that has a finite length change in the absence of an external field.cAFQ L and MSI share some similarities, as both possess a jump discontinuity in the total length change.However, a qualitative distinction is the lack of a jump or peak in quadrupolar contribution for cAFQ L , as compared to MSI.Furthermore, the MSI magnetostriction has a dominant linear-in-h scaling behavior for the MSI before the jump, while the cAFQ L has an overarching nonlinear scaling before S z (0) becomes fully polarized.All of these differences demonstrate the uniqueness of the non-Kramers MSI and MPOs magnetostriction signatures.

V. COMPARISON WITH KRAMERS MAGNETICALLY ORDERED PHASES
A natural comparison is with the more prevalent Kramers ions of the pyrochlore family, which supports the usual mag-FIG.3. Length change, L L , along the (1,1,1) direction under applied [111] magnetic field, h, for the various classically multipolar ordered phases of non-Kramers ions (J z± = 0): (a) coplanar antiferroquadrupolar (cAFQ L ), (b) a second coplanar antiferroquadrupolar (cAFQ L ), and (c) and (d) ferroquadrupolar (FQ L± ) where the ± denote J ±± > 0 and J ±± < 0, respectively.The dashed vertical lines denote regions of discontinuity in the length change, intimately linked to the discontinuity arising from S z (0) becoming fully polarized.The green, blue, and red curves denote the length change arising from the XY pseudospin (quadrupolar), Z pseudospin (dipole), and combined contributions, respectively.Just as in Fig. 2, the dipole plot line is (slightly) purposely shifted from the total length change plot line to more easily visualize the individual contributions.The values of the chosen exchange couplings and lattice-pseudospin couplings are presented in Appendix F. netic dipole moments J x,y,z .Key differences between Kramers and non-Kramers ions are (i) J z± = 0 which causes mixing between the spin-ice and splayed ferromagnet phases and (ii) the magnetic field can couple at linear order to the XY components of the pseudospins for Kramers ions.These differences are a consequence of the pseudospin components being mere dipole moments and are thus odd under time reversal.We present the classical T = 0 phase diagram and the magnetostriction behaviors of the magnetically ordered phases in Appendixes K and L. The classical phase diagram for Kramers ions possesses a variety of broken-symmetry phases: a generalized splayed ferromagnet (SFM), a coplanar antiferromagnetic Palmer-Chalker (PC) phase, and a 1D manifold of antiferromagnetic states.The magnetostriction behaviors for Kramers ions possess jump discontinuities in the length change for increasing field strengths, and for certain phases there are multiple such discontinuities.
We can draw specific contrasts between the associated Kramers and non-Kramers states.For instance, the 1D manifold-like states in Kramers case have a discontinuity, while the corresponding FQ L states of non-Kramers ions undergo smooth length change under increasing field.Analogously, SFM and cAFQ L can be distinguished as SFM has a vanishing length change at zero magnetic field, while for cAFQ L it is finite.Finally, PC and cAFQ L can be differentiated as PC has two discontinuous points in the length change, while cAFQ L has only one.Such key differences in the length change behaviors of Kramers and non-Kramers ions highlight the broad applicability of magnetostriction in pyrochlore materials.Furthermore, since Kramers ions are more commonly examined with conventional probes of magnetic ordering (most notably neutron scattering), magnetostriction can thus serve as useful corroborating evidence.We contrast this with the non-Kramers situation, where there is a dearth of probes available, and where each of MPOs possess distinct features (Fig. 3) that allows each of them to be individually identified (and distinguished from non-Kramers MSI).This comparison thus also serves to emphasize the suitability of magnetostriction to non-Kramers ions.

VI. DISCUSSION
In this work, we proposed that magnetostriction is an ideal probe of multipolar quantum spin-ice ground states and multipolar ordered states in the pyrochlore oxide family.Employing a symmetry-based approach, we constructed an elastic strain coupling to the local rare-earth-metal moments, in the presence of an external magnetic field.We studied all the possible classically ordered states of the pseudospin-1/2 model, and we also employed a 32-site exact diagonalization of the quantum model to examine the multipolar quantum spin ice in non-Kramers ions.We found that the multipolar quantum spin-ice phase in non-Kramers compounds has a unique magnetostriction behavior that is distinct from the other multipolar ordered phases (which themselves have unique behaviors).We note that while we discovered the KI phase under a finite magnetic field for the parent multipolar spin-ice phase, quantum Monte Carlo studies on the simplified X X Z model suggest that the so-called "monopole supersolid" ordered phase could also develop into a KI when subjected to a magnetic field [43].This KI is, however, preceded by the rapid demise of the XY orderings, which may provide an additional feature to the magnetostriction.
Our concrete theoretical results for magnetostriction under [111] and [110] magnetic fields provide a guide for targeted experimental investigations for MSI.Indeed, recently presented experimental data of magnetostriction in Pr 2 Zr 2 O 7 [35] seems to be qualitatively similar to our results.Experimentally examining the length change for [110] magnetic field would be an important next step in verification of our selection rules of magnetostriction.Our study is also broadly applicable to other multipolar quantum spin candidate materials, Pr 2 Sn 2 O 7 [54][55][56] and Pr 2 Hf 2 O 7 [57][58][59].In terms of future work, it would be interesting to examine finite-temperature length-change behaviors, such as thermal expansion.As well, identifying the nature of the finite-field quantum kagome ice state of the generic pyrochlore model (as well its possible connection to the resonating plaquette state in the simplified X X Z model [29,43]) would be an important study.Such studies would provide insight into the nontrivial fractionalized excitations predicted in quantum spin ice, such as the emergent monopoles and photon.To further place our work in the physically relevant context of Pr 2 Zr 2 O 7 , it would also be an interesting and important future study to examine and incorporate the impacts of disorder [11,60] in the context of magnetostriction.Finally, it would also be intriguing to examine the study of magnetostriction in other frustrated lattices (with different symmetries) which are candidates for QSLs.It would be fascinating to explore whether those systems also possess strong magnetostriction signatures, for both their proposed QSL and/or any nearby ordered phases.
Note added.After completing our work, we became aware of a parallel and independent theoretical work by Subhro Bhattacharjee and Roderich Moessner, concerning quantum spin-ice physics of Pr 2 Zr 2 O 7 [61].We thank them for communicating their results to us.The authors declare no conflict of interest.

APPENDIX A: LOCAL BASES AND PSUEDOSPIN MODEL MATRICES
The pyrochlore lattice is an underlying face-centred cubic (fcc) Bravais lattice with four sublattices per unit cell.We define the following local bases on each sublattice α in Table I.
Within the local bases, we employ the γ i j matrix, which has the following matrix representation: where w = e 2πi/3 .

APPENDIX B: SYMMETRY TRANSFORMATIONS OF PSEUDOSPINS, MAGNETIC FIELD, AND ELASTIC STRAIN UNDER D 3d
The D 3d point group can be generated by the two following elements, which written in an orthonormal basis (R 3 space) are where S − 6 is an improper rotation about the z axis by π/3, and C 21 is a π rotation about the y axis.Using these generators, we can transform the pseudospin-1/2 quantities (on sublattice α) as And finally, the elastic tensor transforms in the usual manner, i.e., ← → → A ← → A T , where A is the symmetry element.

APPENDIX C: RELATING QUANTITIES IN LOCAL AXES TO GLOBAL AXES
We present the transformation of the relevant quantities from the local axes to the global axes.The vector-like quantities such as the magnetic field are transformed using h α = P −1 α h, and the tensor strain is transformed using α = P −1 α P α .Here P α is the change of basis matrix for sublattice α; i.e., its columns contain the basis vectors of the given subalttice as denoted in Table I.For concreteness, we present the local-to-global magnetic field transformations below: We reiterate that h x,y,z are magnetic field components in the global basis (i.e., no sublattice index).

APPENDIX D: GENERALIZED MAGNETOSTRICTION EXPRESSIONS FOR NON-KRAMERS IONS
To determine the magnetostriction (or total length change) expression ultimately requires knowledge of the relative length change in terms of the elastic strain components.In that respect, we appeal to the relative length change expression of (as previously derived in Ref. [62]) where i j ≡ 1 2 ( ∂u i ∂x j + ∂u j ∂x i ) is the standard strain tensor (in the global basis) and ˆ i is the ith component of the unit vector ˆ .As is clear from Eq. (D1), determining i j is essential to find the fractional length change.To do so, the strain tensor in Eqs. ( 3) and ( 4) first needs to be rewritten in the global basis, using the change of basis described in Appendix C.This change of basis ensures that all elastic strain-dependent quantities are written in terms of elastic strain normal modes in the global basis μ,ν,xy,xz,yz .The subsequent total elastic free energy, F elastic = F lattice + F XY,NK + F Z , is then minimized with respect to the elastic strain normal modes to yield extremized elastic strain expressions, presented in Appendix E. Finally, the extremized elastic strains are inserted into Eq.(D1) to yield the generalized magnetostriction expressions along the direction of interest, .

APPENDIX E: EXTREMIZED ELASTIC STRAIN TENSOR EXPRESSIONS FOR NON-KRAMERS IONS
Extremizing the elastic free energy with respect to the normal modes ( δF elastic δ i j = 0) yields the expressions in Eqs.(E1)-(E6).We emphasize that the above expressions can be used [in combination with Eq. (D1)] when finding the length change along any direction of interest: In the above Eqs.(E1)-(E6), we use the superscript to denote the extremized elastic strain.We note that the magnetic field has also been rewritten in terms of the global basis as described in Appendix C.

APPENDIX F: NUMERICAL VALUES OF CHOSEN COUPLING CONSTANTS CHOSEN
We take J zz = 1 in this study.For the (multipolar) spin ice (M)SI magnetostriction behaviors, we choose J ± /J zz = 0.02 and J ±± /J zz = 0.05.For the cAFQ L magnetostriction, we choose J ± /J zz = −0.5 and J ±± /J zz = −0.5.For the cAFQ L magnetostriction, we choose J ± /J zz = −0.5 and J ±± /J zz = 0.5.For the FQ L+ magnetostriction, we choose J ± /J zz = 0.72 and J ±± /J zz = 0.5.For the FQ L− magnetostriction, we choose J ± /J zz = 0.72 and J ±± /J zz = −0.5.We take Finally, we take δ 1 = 0.00075 and δ 2 = −0.000088 to emphasize the perturbative nature of the quadratic-in-h magnetic field coupling.The numerical values for the pseudospin-lattice couplings are taken with comparison to an experimental study of Pr-based heavy fermion compound, PrIr 2 Zn 20 [16].PrIr 2 Zn 20 shares similarities with Pr 2 Zr 2 O 7 in that both their interesting phenomena arise from Pr ions' f 2 electrons.Taking the above coupling constants yields magnetostriction behaviors that are of the same scale as the reported study.Indeed, the physical scale of ( L/L) ∼ 10 −6 for the relative length change is also observed in magnetostriction studies in other f electron heavy fermion compounds [63,64], as well as in pressurized Kitaev materials [65].The actual value (or ratio) of the coupling constants can be determined by employing the proposed length change behaviors in conjunction with experimental measurements.As an example, by subtracting off the leading linear-in-h scaling behavior in the experimental length change measurements for the [111] field and (1,1,1) direction allows the determination of (2k 1 + k 2 ) in Eq. ( 7) and subsequently (k 1 − k 2 ) from the (1,1,0) length change in Eq. ( 8), as we have numerically computed the pseudospin configurations.To subsequently extract out the remaining g 1,2,3,4 couplings, it requires additional length change measurements.For example, one may employ the gradients of the aforementioned two [111] field length changes along with the gradients of the two [110] field length changes of Fig. 9 to obtain the lattice-dipole couplings g 1,2,3,4 .

APPENDIX G: CLASSICAL SOLUTION TO PSEUDOSPIN-1/2 MODEL
The classical ground state of Eq. ( 1) is that of a 4-sublattice q = 0 ordering, as Eq. ( 1) can be written as the decoupled sum over an individual (up or down) tetrahedron.Thus, any classical configuration that minimizes the energy of a single (up or down) tetrahedron automatically minimizes the total Hamiltonian.The subsequent magnetic orderings over the entire pyrochlore lattice is that of magnetic orderings repeated over each tetrahedron.Since the non-Kramers ions involve quadrupolar moments, these orderings are in fact multipolar orderings (rather than magnetic orderings as in the Kramers case).

APPENDIX H: EXACT DIAGONALIZATION (ED)
The ED ground state is obtained by employing the quantum lattice model solver package H [66].A central step in this package is to represent the pseudospin operators in terms of fermionic operators, namely Using this formulation in Eq. ( 1), eigenenergies, eigenstates, one-body Green's function, and two-body Green's function are obtained.The one-body Green's function permits the extraction of the expectation value of the pseudospin operator, i.e., S μ i .The two-body Green's function allows the pseudospin-pseudospin correlation function to be obtained, i.e.. S μ i S ν j .The convergence factor of the Lanczos algorithm is determined by the condition of whether the relative error between the ground-state energy at a given step and that of the previous step is less than 10 −9 .

The 16-and 32-site ED cluster study
We present in Fig. 4 the 16-and 32-site ED clusters that we employ in this study, which depict the precise locations of the sublattices with respect to the boundary conditions.The 16-site cluster is composed of two Bravais lattice points in the global x and ŷ direction and a single lattice point in the global ẑ.The 32-site cluster is formed by having two Bravais lattice points in each of the x, ŷ, and ẑ directions.
The computationally less intensive 16-site ED study provides the h = 0 phase diagram in Fig. 5.The location of the ED phase boundaries in the 16-site study guides us in choosing the J ± = 0.02J zz and J ±± = 0.05J zz parameter choice to investigate the 32-site quantum spin ice under a magnetic field.As seen, we confirm the existence of three distinct phases separated by phase boundaries.The phase boundaries are characterized as a singular point in the second deriva-tive of ground-state energy with respect to the two coupling constants, i.e., singularity in ∂ 2 E /∂J 2 ± and ∂ 2 E /∂J 2 ±± .We compare the qualitative similarity of Fig. 5 to the classical phase diagram of Fig. 1 as well as the gMFT phase diagram in Ref. [31].We note that the precise location of the phase boundaries is different when comparing ED to gMFT FIG. 6.The intensity distribution of S 22 (q) for each phases.The parameter sets are chosen as (J ± , J ±± ) = (0.03, 0.1) for phase I, (J ± , J ±± ) = (0.4,0.1) for phase II, and (J ± , J ±± ) = (0.1, 0.8) for phase III, respectively.
(or even classical) studies.Indeed, the ED phase boundary location along the J ±± = 0 axis is at J ± = 0.056J zz .This phase boundary is much closer (compared to gMFT) to the quantum Monte Carlo phase boundary of the X X Z model of J ± = 0.052J zz [67].This comparison provides strong support to the validity of our ED findings.Moreover, the 32-site ED favorably compares to other numerical techniques such as numerical linked cluster methods in pyrochlore material Yb 2 Ti 2 O 7 [68].
To understand the nature of these phases, we examine the static-spin structure factor.The intensity of the structure factor is defined as where μ sums over the three components {x, y, z} of the pseudospin, α and β are sublattice indices {0, 1, 2, 3}, N s is the total number of sites, and i, j are site locations of sublattice α, β, respectively; in the N s = 16 (N s = 32) site cluster, there are four (eight) such i, j locations each.The wave number k is represented by using primitive reciprocal vectors b i as k = 3 i=1 q i b i .From this notation, we can easily notice that the first Brillouin zone is for −1/2 < q i < 1/2.We explicitly note that for the 16-site cluster, there are two momenta points in the k x and k y direction, namely that of 0 and π , while the k z direction only has one momentum point of 0, since there is only one Bravais lattice point in the z direction.
We present in Fig. 6 the static structure factor S 22 (q) for each of the three regions of Fig. 5, which provides information on the long-range correlation effects in the cluster.In particular, the location of the peak structure provides information as to the nature of the multipolar order realized in the system.For region I, we find intensity peaks at q = (0, 0) and q = (π, π ), which is a reflection of a lack of an ordering wave vector within the 16-site cluster.This lack of order seems qualitatively consistent with a QSL phase.On the other hand, phases II and III have a single peak located at q = (0, 0), which validates a q = 0 ordered state.
To distinguish phases II and III is, however, challenging as the expectation value of the local pseudospin moment in the absence of a symmetry-breaking magnetic field is always zero in ED.By studying the pseudospin-pseudospin correlation function, we can fortunately demonstrate the consistency of these phases with the classical ordered FQ L and cAFQ L phases.For instance, in phase II, the nearest-neighbor correlation is S x(y) α S x(y) β > 0 and S z α S z β > 0 + .Here 0 + indicates a positive number that is an order of magnetitude smaller than S x(y) α S x(y) β .This indicates a ferrolike correlation in the xy local moments (that is dominant over any z-component correlation), which is consistent with the FQ L .Similarly, for phase III, the nearest-neighbor correlation is consistent with the cAFQ L phase.Thus, we can reasonably claim that the 16-site ED phase diagram matches well with the expected phase diagram from gMFT.We subsequently use the 16-site ED phase diagram as a guide for the choice of parameters to use for the 32-site ED investigation of quantum spin ice.In particular, in anticipation that the phase boundaries will likely shift, we choose J ± and J ±± to be deep in phase I (the likely quantum spin ice phase) and away from the phase boundaries.

Extracting ED pseudospin expectation values
An inherent challenge in ED studies is in extracting the pseudospin expectation values, as spontaneous symmetry breaking is only captured in the thermodynamic limit [69].As such, in finite-sized clusters (in zero magnetic field), the pseudospin expectation values are always zero.This issue can be avoided in a magnetic field, as the coupling of the field to the dipole moment explicitly breaks the symmetry, thus rendering a finite dipole expectation value.However, in the absence of a field, we sketch the general strategy that can be used, where we employ the two-point pseudospin correlators and the intensity of the structure factor to divine the pseudospin expectation values.In particular, we use the numerically produced S z i values as a benchmark to extract out the pseudospin expectation values of the x and y components.Even though the quadrupolar contribution to the length change only involves the X and Y pseudospin components on sublattices 1, 2, and 3 (as will become clear below), it is helpful to have information of the ordering on sublattice 0.
We first recall that the pseudospin on sublattice 0 is fully polarized under a [111] field to yield a trivial ordering on the (as discussed in the main text) yields an averaged (over all degenerate states) local pseudospin configuration.

APPENDIX J: CLASSICAL NON-KRAMERS SPIN ICE IN [110] MAGNETIC FIELD
We present in Fig. 9 the classical magnetostriction behavior along the (1,1,1) and (1,1,0) directions for a [110] magnetic field.As seen, there is a lack of any clear or distinct features.The reason lies with the fact that both the magnetic field and the ice rules can be simultaneously satisfied for this field direction.For the [110] field, the magnetic field couples solely to sublattice 0, 3 and fails to do so to sublattices 1, 2; i.e., only ĥ • (x, ŷ, ẑ) 0,3 = 0, while ĥ • (x, ŷ, ẑ) 1,2 = 0.As such, pseudospins on sublattices 0 and 3 respectively get aligned parallel (+ẑ 0 ) and antiparallel (−ẑ 3 ) to the field, while the other sublattices conspire together to satisfy the ice rules, i.e., partial degeneracy of the ice rules remains, with sublattices 1, 2, taking S z (1) = { 1 2 , − 1 2 } and S z (2) = {− 1 2 , 1 2 }, respectively.This situation is valid for any finite field values.Thus, there is no transition to any fully polarized state in the large field limit and no observable transition (unlike the [111] direction).

APPENDIX K: CLASSICAL PHASE DIAGRAM OF KRAMERS IONS
The classical Kramers phase diagram of Fig. 10 provides a variety of possible phases: a blended phase composed of spin ice (SI) and splayed ferromagnet (SFM) of the same T 1g symmetry, a coplanar antiferromagnetic Palmer-Chalker (PC) phase of T 2g symmetry, and a 1D manifold of states with E g symmetry.An obvious distinction between the Kramers and non-Kramers phase diagram is that SI and SFM phases blend together for Kramers ions, while the corresponding non-Kramers phases are separated by a phase boundary for non-Kramers ions.This is a consequence of J z± = 0, which allows the two aforementioned order parameters to mix.One can easily notice this by expressing Eq. ( 1) in terms of classical order parameters (orderings) on a single tetrahedron, where we use the definition of the order parameters as presented in Appendix N, and we drop the gerade subscript for the order parameters for brevity.From the last term in Eq. (K1), the two T 1 symmetry magnetic orderings mix with each other when J z± = 0. Consequently, for Kramers ions, there exists a common region in the phase diagram with coexisting SI and SFM ordering.Depending on the location in this common region, the state is more SI-like or more SFM-like, which we label as SI dominant and SFM dominant, respectively.

APPENDIX L: MAGNETOSTRICTION EXPRESSIONS OF KRAMERS PYROCHLORE MATERIALS
We now turn to examining the magnetostriction behavior of classically ordered Kramers phases.Since all the pseudospin components are magnetic dipole moments, they all couple to an external magnetic field at linear order, H mag,K = −h •

)FIG. 4 .
FIG. 4. The 16-and 32-site ED clusters.[(a), (c)] Schematic of the respective 16-and 32-site unit cells with the yellow (blue) planes denoting the kagome (triangular) planes.[(b), (d)] Top view of the respective 16-and 32-site clusters from the [111] direction.The four different colors denote sites on the four different triangular and kagome layers in panels (a) and (c), respectively.The periodic directions are denoted by the axes.

FIG. 8 .
FIG. 8. Local pseudospin configuration for non-Kramers SI in a [111] magnetic field.The degenerate kagome spin ice phase is denoted by the yellow shaded region, stable quantum spin ice is indicated by the orange shaded region, and the fully polarized state is depicted by the indigo-blue shaded region.Solid lines (unfilled squares) indicate classically (ED) obtained local pseudospin configurations on each sublattice.

FIG. 11 .
FIG. 11.Length change, LL , along the (1,1,1) direction under applied [111] magnetic field, h, for the various classically magnetically ordered phases of Kramers ions (J z± = 0.25J zz ).The dashed vertical lines denote regions of discontinuity in the length change, intimately linked to the discontinuity in the pseudospin expectation values.The green, blue, and red curves denote the length change arising from the XY pseudospin (quadrupolar), Z pseudospin (dipole), and combined contributions, respectively.
Center for Quantum Materials at the University of Toronto.Y.B.K. is supported by the Killam Research Fellowship of the Canada Council for the Arts.Computations were performed on the Niagara supercomputer at the SciNet HPC Consortium and on the Graham supercomputer of Compute Canada.SciNet is funded by the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund, Research Excellence; and the University of Toronto.We thank the Center for Advanced Computation at the Korea Institute for Advanced Study for providing computing resources for this work.S.B.L. is supported by the National Research Foundation Grant No. NRF2017R1A2B4008097. M.H. is supported by the Japan Society for the Promotion of Science through the Program for Leading Graduate Schools (MERIT) and Overseas Challenge Program for Young Researchers.We thank Satoru Nakatsuji, Mingxuan Fu, and Nan Tang for helpful discussions on their experimental data on Pr 2 Zr 2 O 7 and sharing their insights.Y.B.K. conceived and supervised the research.A.S.P. and M.H. performed the calculations in this work.S.B.L. provided important preliminary results regarding the pseudospin-lattice couplings and contributed to the discussion of the overall results.All authors contributed to writing the paper.

TABLE I .
Local sublattice basis vectors.
[36]e we include the nonvanishing g-tensor components g xy , g zz .As a typical example, we take the estimated g-tensor values of Yb 2 Ti 2 O 7[36]: (g xy , g zz ) = (4.18,1.77).Due to the magnetic dipole nature of the XY moments, we have in addition to Eq. (4),