Disorder and Quantum spin ice

We report on diffuse neutron scattering experiments providing evidence for the presence of random strains in the quantum spin ice candidate Pr2Zr2O7. Since Pr is a non-Kramers ion, the strain deeply modifies the picture of Ising magnetic moments governing the low temperature properties of this material. It is shown that the derived strain distribution accounts for the temperature dependence of the specific heat and of the spin excitation spectra. Taking advantage of mean field and spin dynamics simulations, we argue that the randomness in Pr2Zr2O7, promotes a new state of matter, which is disordered, yet characterized by short range antiferroquadrupolar correlations, and from which emerge spin-ice like excitations. This study thus opens an original research route in the field of quantum spin ice.

In condensed matter physics, disorder usually tends to freeze the degrees of freedom. This is for instance the case in the Anderson localization phenomenon [1] or in the spin glass transition [2]. In frustrated systems, disorder is also expected to dramatically impact the ground state, possibly engendering new states of matter. For instance, in pyrochlore antiferromagnets, considered as the archetype of 3D geometrically frustrated systems, weak structural and magnetic disorder can induce spin-glass [3,4] or even "topological spin glass" phases [5] in spin ice. When the magnetic moment in the latter materials is formed by non-Kramers ions, intrinsically very sensitivity to strain, disorder or local distortions, as for instance in Tb 2 Ti 2 O 7 , Tb 2 Hf 2 O 7 , Pr 2 Zr 2 O 7 or Pr 2 Sn 2 O 7 , the physics is even more intriguing. A bunch of experimental studies reported the lack of long range order down to very low temperature along with an unexpected spin fluctuation spectrum, indicating a strong density of low energy excitations [6][7][8][9][10][11][12][13]. It was then suggested that instead of driving glassy behavior, disorder may open a new route in stabilizing the enigmatic "quantum spin ice" state [12,14].
While this scenario appears extremely appealing, raising eventually the hope of tuning new exotic states of matter by disorder, without even touching the magnetic lattice, it remains to show that the actual disorder in real materials is indeed consistent with these ideas. In this work, we address these issues in the particular case of Pr 2 Zr 2 O 7 which is ideally documented among quantum spin ice candidates. We first report on neutron scattering measurements providing evidence for an intrinsic distribution of strains, ensuring a direct experimental foundation for our point. We then examine in details the consequences of this randomness on the low temperature The strain induced by a defect in the pyrochlore structure is sketched using red concentric spheres. The CEF scheme [11,15] is shown schematically on the left side, for the unperturbed (blue) and perturbed (red) sites. Double (resp. single) blue lines correspond to doublets (singlets). The shaded red rectangles illustrate the broadening of the CEF level by disorder. Finally, the black lines indicate the local CEF axis directions. magnetic properties. It is found that this random strain distribution indeed accounts for a number of experimental features. Eventually, we find that it corresponds to a strong perturbation. Pr 2 Zr 2 O 7 then hosts a competition between disorder and the "native" effective interactions, namely the Ising coupling J zz (considered in Ref 14) and the transverse quadrupolar terms J ± , known to play a prominent role in Pr 2 Zr 2 O 7 [15]. With the support of real space mean field and time-evolving spin dynamics simulations, we show that the competition with J ± especially promotes a new kind of spin liquid, disordered (or weakly ordered), characterized by short-range antiferroquadrupolar correlations and from which emerge a peculiar spin-ice like excitation spectrum. This study takes advantage of an abundant literature on Pr 2 Zr 2 O 7 . As sketched in Figure 1, the physics is governed by the ↑ and ↓ degenerate states of the rare earth crystal field (CEF) ground state doublet. The corresponding magnetic moments sit on the vertices of the pyrochlore lattice (made of corner sharing tetrahedra) and point along local 111 directions (see Figure 1). In Pr 2 Zr 2 O 7 , the {|↑ , |↓ } subspace is well protected since the first CEF excited state is located at 10 meV [11,13,16,17]. Furthermore, the Curie-Weiss temperature inferred from magnetic susceptibility is negative, indicating antiferromagnetic interactions. The specific heat shows a broad peak at about 2 K [11,13,16,[18][19][20], bearing some resemblance with the classical spin ice Dy 2 Ti 2 O 7 , along with an upturn at low temperature attributed to the hyperfine contribution. Above 1 T, the broad anomaly shifts to larger temperatures [15]. Finally, the magnetic excitation spectrum S(Q, ω), measured by inelastic neutron scattering, has been described as a dynamical spin ice mode consisting of a broad inelastic response centered around ∆ ≈ 0.4 meV. Its structure factor strongly resembles the famous spin-ice pattern, with pinch points and arm-like features along hhh directions [11,12,15].
In the particular case of non-Kramers ions, the degeneracy of the ground doublet is easily lifted by perturbations, defects, local deformations and strains which, by virtue of the magneto-elastic interaction, perturb the electronic density over the rare earth sites. The Ising |↑↓ doublet turns into "tunnel-like" wave-functions [8][9][10]12] and split by a random quantity ∆ which directly reflects the strength of the perturbations (see Figure 1). In this local picture, as we show below, the specific heat consists in multiple Schottky anomalies due to the distribution of splittings and the large width of the spin excitations spectrum arises from the transitions within the randomly split doublets. These transitions can be observed by neutrons because their cross section, proportional to a|J z |s = ↑ |J z | ↑ , is non-zero. This is at variance with the unperturbed case since the non-Kramers character imposes | ↑ |J| ↓ | ≡ 0.

Lattice strain and diffuse scattering
Pr 2 Zr 2 O 7 samples have been prepared as carefully as possible, and especially annealed to prevent the presence of Pr 4+ ions, resulting in a beautiful green transparent color [21]. Their crystalline structure is very well refined by a perfect pyrochlore structure, the amount of defects being beyond the accuracy of standard diffraction, typically 1% [15,17]. The pyrochlore structure of Pr 2 Zr 2 O 7 actually contains two types of cations, Pr 3+ and Zr 4+ on 16d and 16c sites respectively. The two oxygen sites O1 (48f) and O2 (8b) are fully occupied, whereas the 8a site is vacant.
Lattice strains in Pr 2 Zr 2 O 7 are first evidenced by means of polarized neutron scattering measurements. Fig. 2 shows Non-Spin-Flip (NSF) diffuse scattering maps recorded in the (HH0, 00L) scattering plane (see Supplemental Information). Anomalies appear in the vicinity of the nuclear Bragg reflections (222), (311) and (400), taking the form of butterfly features elongated along specific directions (Fig. 2). In contrast, the nuclear Bragg peaks (111) and (220) show quasi isotropic shape with negligible diffuse scattering. The measurements were repeated at three temperatures (10 K, 1 K and 50 mK). No significant evolution could be noticed, hence suggesting that this diffuse scattering is associated with structural disorder and lattice strains.
In many ordered pyrochlores, the less energetic defects mechanisms [22] are: (i) cation antisite disorder, namely substitution of R 3+ cations by T 4+ ones or the opposite; (ii) creation of anion Frenkel pair, consisting of a vacant 48f site and an interstitial oxygen occupying a 8a site. The association of both processes is often at play. It influences electronic or ionic conductivity, and helps the accommodation of radiation damages and the resistance to amorphization [23]. However, these disordered fluoritelike defects are not the only source of lattice strains in the ordered pyrochlore structure. Other mechanisms involving lattice displacements could be at play. For example, tiny or short range ordered distortions lead to dimerlike magnetic excitations in GeCo 2 O 4 ordered spinel [24] as well as a kind of frustrated ferroelectricity in niobium based pyrochlores [25]. In Pr 2 Zr 2 O 7 , a local off-centering of the Pr 3+ ions from the ideal pyrochlore 16d sites has been reported [26] as the dominant mechanism, similar to observations for the trivalent cation in Bi 2 Ti 2 O 7 [27] and La 2 Zr 2 O 7 [28].
Anyhow, irrespective of the actual nature of disorder, the lattice strains can be analyzed using the same formalism. The diffuse scattering maps have been worked out in the framework of the so-called Huang scattering [29] (see Supplemental Information), assuming random static displacements of the atoms of the structure. The scattered intensity at a reciprocal lattice vector Q = τ + q around a Bragg peak τ is written: where I bg is a constant background, |F (τ )| the structure factor of the Bragg reflection, S Bragg (τ ) a normalized two-dimensional Gaussian and c is the concentration of defects leading to diffuse scattering modeled by the func-  tion S Diff (Q). Introducing the Fourier transform s (q) of the displacement field s (r), the diffuse scattering can be modeled by the following function: (3) which is composed of symmetric |τ · s (q)| 2 and antisymmetric iτ · s (q) parts and involves the static Debye factor of the defects L τ (q). In the following, we assume that the strain field s (r) is isotropic around a defect (see [30,31] along with Supplemental Information). The disorder is parametrized by spherical defects of typical volume V C creating an internal local pressure on the lattice P 0 . We assume V C = 10.92Å 3 , corresponding to the difference in atomic volume of Pr 3+ and Zr 4+ ions and P 0 = 10 eV, which is a typical value in a number of materials [32,33].
Around each Bragg peak, the diffuse scattering was described by the model given by Eq. (2) using the above fixed parameters V C and P 0 , and taking the structure factor F (τ ) and the concentration of defects c as fitting parameters. Note that in the course of our analysis, we neglect the effects of sample mosaic spread and imperfect instrumental resolution. Indeed, a two-dimensional Gaussian fit of the nuclear Bragg peaks yields widths much smaller than data pixelation, rendering a convolution of the calculated I (Q) with the instrumental reso-lution function unnecessary (see Supplemental Information).
Examples of fits are given in Fig. 2 for selected Bragg peaks, together with the calculated part of the Huang diffuse scattering. The fitted structure factor of the Bragg peaks nicely corresponds to calculations of the perfectly ordered structure. The concentration of defects can be estimated to c = 1.1(5) · 10 −3 , when averaged over the 7 measured Bragg peaks. This yields an average volume change induced by one defect of ∆V /V 1.3 · 10 −3 , namely an order of magnitude hardly measurable by conventional diffraction (but which could be accessed by neutron Larmor diffraction, which allows measuring the intrinsic width of the unit cell volume distribution [34,35]).
According to Refs. 36 and 37, the above model corresponds to a random strain field described by the distribution function: where e = {e 1 , e 2 , e 3 , e 4 } denotes the four independent components of the deformation tensor, e 2 = i=1,..,4 e 2 i is the total strain and γ is a material dependent dimensionless width that reflects the level of defects. Following Malkin et al. [37], γ can be computed from the unit cell volume V dependence on point defect concentration c using the relationship: where σ is the Poisson ratio. Unfortunately, V (c) is unknown for Pr 2 Zr 2 O 7 so that γ cannot be accurately determined from the above neutron data only. However, an analysis of existing data on the sister compound Tb 2 Ti 2 O 7 [38] (see Supplemental Information) yields an estimate of the order of magnitude of γ, namely 10 −5 −10 −4 . Inputing the above refined values for ∆V /V and c for our Pr 2 Zr 2 O 7 sample in Eq. 5, we indeed get γ 3.2 ± 1.5 · 10 −5 . This polarized neutron study thus highlights the existence of strains in the material, not detectable with standard diffraction measurements. With this result in hand, we now examine the consequences of those random strains on the low temperature magnetic properties.

The picture of isolated doublets
As a first step towards this objective, we consider a simplified picture, where the Pr 3+ are supposed independent, yet split by a random perturbation due to the strain. Using, at each site, a pseudo spin 1/2 (σ x , σ y , σ z ) spanning the {|↑ , |↓ } subspace, the magneto-elastic interaction takes the simple following form (see Supplemental Information): where v i is an auxiliary random field that depends on the strain distribution g(e). The real and imaginary parts of the random field v are characterized by an average v = 0 and a standard deviation (which physically corresponds to the strength of the disorder) δv ≈ kγ (7) where k is a numeric constant which merges the magnetoelastic coefficients (k ∼ 1.6 × 10 4 K, see Supplemental Information). At this simplified level of approximation, the model given by Hamiltonian (6) describes an assembly of uncoupled two-levels subsystems. The | ↑ and | ↓ states recombine in |a and |s , split by the random variable ∆ = 2|v|. The shape of the probability density p[∆] is shown in Fig. 3a. It exhibits an asymmetric profile with a maximum at ∆ m ≈ 4/3 δv ≈ 4/3 kγ.
In the two following paragraphs, the specific heat along with S(Q, ω), calculated on the basis of Hamiltonian (6) and weighted according to p [∆], are compared to exper-iments.

Specific heat
The specific heat C p was measured in Ref. 11 from very low temperature up to 25 K as well as in Ref. 15 as a function of field. It shows a broad peak between 0.3 K and 10 K, followed at very low temperature by a steep upturn attributed to hyperfine effects due to the nuclear Pr 3+ spin.
In the above model of independent split doublets, the specific heat is written: Our simulation of the specific heat C p (T ) is shown in Fig. 3b (black solid line, together with the 1/T plot of ln C p (T ) in inset). It consists of multiple Schottky peaks due to the superposition of individual random splittings ∆ induced by local strains. This contrasts with the interpretation suggested in Ref. 11 in terms of monopole quantum dynamics. We find that for a value γ = 1.25 ×10 −4 (hence within the uncertainty range stimated from Huang scattering), corresponding to a disorder strength δv ≈ 2 K, C p (T ) shows a peak at 2 K, matching the experimental result, with a ( 25 %) overestimate of the peak amplitude. The entropy released by the calculated specific heat anomaly is expectedly Rln 2, which is larger than the value estimated from the 4f contribution extracted from the data [11], and in agreement with the higher calculated peak value.
In Ref. 11, it was proposed that the upturn observed at lower temperature could be tentatively attributed to a hyperfine Schottky anomaly. However, this work pointed out that, provided the ground Pr 3+ doublet remains degenerate, this hypothesis would lead to a too large upturn, the specific heat showing then a maximum at 0.14 K peaking at 7 JK −1 (mol.Pr) −1 . We show here that the disorder not only allows to qualitatively explain the peak at 2 K, but also resolves this discrepancy at very low temperature. The simulation is performed by adding the hyperfine contribution A I z J z to the Hamiltonian, where A is the magnetic hyperfine constant A 0.055 K [39] and I z is the Pr nuclear spin with quantum number m = 5/2, 3/2, . . . , −5/2. The calculated specific heat is represented as a red curve in Fig. 3b, revealing that the hyperfine upturn has the correct order of magnitude. It is worth mentioning that the magnitude of this upturn directly results from the distribution of splittings: when the strain is small, the ground Pr 3+ doublet survives and remains degenerate, with an Ising behavior. This gives rise to a large hyperfine anomaly of the specific heat. In contrast, when the strain is large, the ground state is a singlet and the hyperfine anomaly weakens. The existence of a splitting distribution thus creates an intermediate situation between these two limiting cases, and finally explains why the upturn of the specific heat at low temperature remains moderate. In this sense, the fact that the upturn remains somewhat underestimated indicates that δv ≈ 2 K is likely an upper limit of the actual strain distribution. We also notice that the amplitude of the low temperature up-turn strongly depends on the sample synthesis, namely the growth rate of the crystal [26].

Spin excitations and Inelastic neutron scattering
Low energy inelastic neutron scattering spectra taken from Refs. 15 and 40, are represented in Fig. 3c at selected temperatures and for a momentum transfer Q=(1.8,1.8,0). They are part of the dynamical spin ice mode reported initially in Ref. 11 and consist of a broad inelastic line whose width increases and whose peak intensity decreases upon heating. Both quantities saturate above 10 K up to 50 K.
In the model of independent doublets defined by Eq.
where 1 + n(ω) = [1 − exp(−ω/k B T )] −1 is the detailed balance factor and G(∆, Γ) represents energy profiles characterized by differences of Lorentzian line-shapes centered at ±∆ with a half width at half maximum (HWHM) Γ. S(Q, ω) thus consists of a series of modes at the random energies ∆ and the global width of the signal is then a direct signature of the randomness. ζ is defined as | a|J|s | 2 and owing to the actual CEF, ζ = 3.4. The neutron cross section is finally given by where a is an instrumental scaling factor, bg(ω) a linear background term and R(Q, ω) the TAS resolution function (see Supplemental Information). It should be stressed here that the non zero cross section is due to the recombination of the ↑ and ↓ states. Otherwise, it would be zero because of the non-Kramers nature of the Pr 3+ ion which imposes ↑ |J| ↓ ≡ 0.
The calculated cross sections at different temperatures are represented in Fig. 3c with the same p[∆] as the one used to calculate the specific heat. As expected from an ensemble of non-degenerate doublets, the overall scattering is inelastic, with a peak around 0.4 meV at the lowest temperatures. Comparing to the experimental data, the spectra at 1.5, 3, 5, 10 and 20 K are indeed satisfactorily reproduced by a simultaneous fit of Eq. 10 to the whole set of curves. Since the physics probed by inelastic neutron scattering is assumed to be driven by random strains, hence by p[∆], it is actually natural for a to be T −independent. This agreement implies that the model captures the physics at play. In turn, the thermal increase of the line width is thus due to scattering by a distribution of electronic transitions with energies being of the same magnitude as T .
The above results advocate that the picture of independent doublets split by the actual disorder can account for a number of features not understood so far: provided that the strength of the disorder δv is of the order of 1 K (with a typical strain γ ≈ 1.25 × 10 −4 ), the temperature dependence of the specific heat is captured as well as the broadening of S(Q, ω). Note that this value of δv should be considered as an estimate only: indeed, on the one hand, γ likely depends on the sample preparation; on the other hand, δv also depends on the exact values of the magneto-elastic coefficients via the coefficient k (see Supplemental information). Furthermore, the exact profile of p[∆] may be more complex. We note for instance that the distribution of splitting proposed by Wen et al [12] on the basis of INS data, is similar to p[∆], yet reaches a maximum at ∆ = 0. It would be interesting to see how this affects C p for instance.
It remains that as far as the doublets are independent, no particular feature is expected in reciprocal space, so that S(Q, ω) should not be Q dependent in this approach. Experimentally, though, S(Q, ω) displays an inelastic spin-ice like pattern and is thus Q-dependent. A thorough description of the actual ground state of Pr 2 Zr 2 O 7 should then incorporate both the randomness and interactions between the pseudo-spins.

The picture of coupled doublets
To go in this direction, we follow Ref. [41][42][43] to take couplings into account. The magnetic and quadrupolar degrees of freedom are described in a unified framework by considering the components (σ x , σ y , σ z ) of the pseudo spins 1/2 that span the {|↑ , |↓ } subspace. Provided these states are well protected from the first excited state, the Ising magnetic moments J z,i , pointing along the 111 z i axes, as well as the quadrupoles can be mapped onto the (σ x , σ y , σ z ) Pauli matrices (see Supplemental Information) via the q ,⊥ coefficients: Q yz = (J y J z + J z J y )/2 = q ⊥ σ y g J = 4/5, g is an effective anisotropic g factor, and, together with q ,⊥ , depend on the actual CEF scheme (g /g J ≈ 6.8, ζ ≈ 3.4, q ⊥ ≈ 0.055 and q ≈ 2.17). On this basis, a bilinear Hamiltonian has been proposed [41][42][43][44][45][46][47]: The γ ij parameters are unimodular matrix coefficients defined in Ref. 44. The upper and lower panels display the magnetic and quadrupolar correlation functions C zz and C ± as a function of J ±,z and for different disorder strength δv. Without disorder (δv = 0), these quantities provide a picture of the mean field phase diagram. This phase diagram encompasses two dipolar phases, spin ice (SI) and "all−in -all−out" (AIAO) as well as two quadrupolar ordered phases, Q-SI and Q-AIAO. In the latter, the pseudo spin is perpendicular to the zi axis. The rectangle displays the region of the phase diagram proposed for Pr2Zr2O7 in Ref. 15. According to our definitions, the σ z are of the same sign in the AIAO states and staggered in the SI phase. C zz then reaches its maximum in the long range ordered phase with C zz = 1/2 × 6 × 1/2 = 3/2 and C zz = 1/2 × (−4 + 2) × 1/2 = −1/2 for the AIAO and SI phases respectively. Similarly the σ x have the same sign in the Q-AIAO states and are staggered in the Q-SI phase. Again, the maximum value is reached in the long range ordered phase with C ± = 3/2 and C ± = −1/2 respectively. The gray rectangles show the range of parameters proposed for Pr2Zr2O7. The red rectangle emphasizes the results for the value of δv estimated from specific heat and neuton scattering.
is indeed simplistic. The full Hamiltonian H+H m−el (Eq (11) and Eq. (6)) then describes an anisotropic Heisenberg model in a random transverse field v i whose distribution is ruled by γ.
The exact solution to this problem is beyond the scope of the present work. We shall however carry on with a simplified picture and consider the mean field approximation (assuming J ±± = 0 as in Ref [15]): At this level of approximation, the model describes an assembly of coupled two-levels subsystems. The eigenstates |± i (at a given site i) are split by the (positive) random quantity ∆ i = 2 |V i | 2 + η 2 i /4 and can be written in a convenient way using local spherical angles φ i and θ i : The |± i states can be seen as intermediate states between the "tunnel-like" wave-functions |s, a given by Eq.
1 and the original | ↑↓ doublet. φ i and θ i have a natural physical interpretation: ± 1 2 cos θ i is the dipolar magnetic moment at site i and ± 1 2 sin θ i (cos φ i , sin φ i ) is the "quadrupolar" moment at the same site. Their values as well as whether |+ or |− is the ground state depend on the parameters J ± , J zz and on the v i random field.
In the absence of disorder, the mean field phase diagram consists of long range ordered magnetic and quadrupolar phases (see Ref 15 along with the left column of Figure 4). It encompasses an antiferromagnetic "all−in -all−out" phase (AIAO) for J zz ≤ 0 (the four spins of a tetrahedron pointing out or towards the center), an ordered spin ice phase (SI) for J zz ≥ 0 (two spins pointing out and two spins towards the center), and two quadrupolar phases denoted Q-SI for J ± ≤ 0 and Q-AIAO for J ± ≥ 0. In the magnetic phases, the pseudo spin are aligned along the local z i axes (θ = 0, π). In contrast, in the quadrupolar phases, the pseudo spin 1/2 are tilted away and order within the XY plane (θ = π/2).
On top of this fully ordered background, the pseudo spins tend to get additional random σ x and σ y components owing to strain. Clearly, the net ordered moment along the z i axes or within the XY plane will decrease with increasing disorder and may even be averaged to zero for strong disorder. Depending on the values of J zz and J ± , magnetic (dipolar) or quadrupolar correlations may or may not survive. It is then convenient to discuss the phase diagram using the magnetic and quadrupolar correlations functions: the sum over " j i " corresponds to a sum over the neighbors connected by J zz or J ± . Figure 4 shows C zz and C ± as a function of J ± and J zz for different disorder strength. Positive (resp. negative) values of C zz indicate AIAO (resp. SI) correlations. In a similar way, positive (resp. negative) values of C ± indicate Q-AIAO (resp. Q-SI) correlations. With increasing disorder, the locations of the different long range ordered phases are still visible but all the more smeared out that J zz δv or J ± δv. In the latter regions of the phase diagram, the long range ordered phases are unstable, being replaced by disordered phases where short range correlations survive (note that the AIAO correlations have a better resistance to disorder than the SI ones essentially because the molecular field is larger in the former).
For δv = 2 K and the coupling parameters proposed for Pr 2 Zr 2 O 7 (J ± = 0.7 K and J zz = −0.5 K), the Q-AIAO order parameter is already suppressed (see Figure 5a), but antiferroquadrupolar correlations on neighboring sites are found to survive (see Figure 4). However, owing to uncertainties on the precise value of δv, we should not confine ourselves to δv = 2 K: from the phase diagram considered as a whole, it can be inferred that Pr 2 Zr 2 O 7 locates in a region where the ground state is disordered or just about to order in the long range Q-AIAO state. Calculations of S(Q, ω) based on real space time evolving spin dynamics simulations (see methods) [48] have been carried out for different δv. For δv = 0, the calculations reproduce the dynamical spin-ice mode observed at 0.4 meV initially reported in Ref. 11, with its typical arm-like features in reciprocal space. The spectra along (h, h, 2) shown in Figure 5b illustrate that with increasing disorder, the energy width of the response grows rapidly, resulting in spectra similar to the ones shown in Figure 3c once convolved with the experimental resolution. The dispersive features also become difficult to distinguish. By integrating S(Q, ω) within 0.38 ≤ ω ≤ 0.48 meV, we obtain the Q-dependence of the dynamical spin ice mode show in Figure 5c. The important point is that the arm-like features, although weaker and marked downturn, are still visible for strong disorder, i.e. in the disordered phase. The pinch points are also considerably blurred, which is consistent with the experiments.

DISCUSSION AND SUMMARY
As detailed in this paper, we thus propose that the randomness in non Kramers ions based pyrochlore magnets (such as Pr 3+ and Tb 3+ ), promotes the rise of new disordered phases, characterized by short range magnetic or quadrupolar correlations. More precisely, by virtue of the magneto-elastic coupling, strains split the non-Kramers doublet of the rare earth ion in a random fashion, recombining the ↑ and ↓ Ising wavefunctions in "tunnel-like" states. J ± and J zz compete with those disorder effects, fostering correlations among dipoles and/or among quadrupoles, deep in the disordered phases, and stabilizing magnetic and/or quadrupolar orderings, at the mean field level, for weak disorder.
Our estimation of the disorder strength δv in Pr 2 Zr 2 O 7 , based on the analysis of polarized diffuse scattering maps and on the fit of various experiments, gives a value of the same order of magnitude as J ± and J zz , locating this material in an intermediate to strong coupling regime.
This scenario provides a very likely and qualitative explanation for a number of experimental features reported in Pr 2 Zr 2 O 7 , such as the lack of long range order, the temperature dependence of the specific heat (including the peak at 2K as well as the upturn at lower temperatures), the energy broadening of the spin excitations spectrum observed by neutron scattering and its specific spin-ice like Q-dependence.
It is worth noting that a related model has been proposed in Ref 14, with J ± ≡ 0 though, hence neglecting quadrupolar effects, the disorder being modeled by a transverse random field h, with averageh and standard deviation δh. Based on a sophisticated treatment, this study points out that disorder provokes quantum superpositions of spins throughout the system, entailing the rise of an emergent gauge structure along with a set of fractional excitations. It also predicts the existence of a "Coulombic Griffiths phase" for large δh and h = 0. Interestingly, this case resembles, to some extent, to the J ± = 0 case considered in the present work. Indeed, at the mean field level, in the Q-AIAO phase, V i = v i − j J ± σ − j behaves as the random field h with a non-zeroh average. As a result, even if dedicated theoretical developments are necessary for a rigorous justification, Pr 2 Zr 2 O 7 could be a relevant candidate for this new state of matter called "Coulombic Griffiths phase".
A few experimental facts remain, however, open questions. For instance, the low temperature upturn of the magnetic susceptibility [13,15] as well as the, although extremely weak, spin-ice elastic response [11] are not captured, at least with these coupling parameters. We emphasize, however, that a positive J zz could help understanding these issues, since such a value would lead to correlated magnetic moments in the spin-ice man- ner. Additional experiments are required to shed light on this point. Furthermore, it would be fruitful to consider other forms of disorder, as for instance the bond disorder pointed out in a recent study of Tb 2 Hf 2 O 7 [13].
As it opens an original research route in the field of quantum spin ice and spin liquids, we hope that this work will arouse new theoretical developments taking especially into account the role of J ± . We also anticipate that experiments on other non-Kramers quantum spin ice candidates as Pr 2 Sn 2 O 7 [49], Pr 2 Hf 2 O 7 [50] and Tb 2 Hf 2 O 7 , where a strongly fluctuating Coulomb liquid phase coexists with defect induced frozen magnetic degrees of freedom [13], will be extremely interesting to test these ideas.