Intermediate magnetization state and competing orders in Dy2Ti2O7 and Ho2Ti2O7

Among the frustrated magnetic materials, spin-ice stands out as a particularly interesting system. Residual entropy, freezing and glassiness, Kasteleyn transitions and fractionalization of excitations in three dimensions all stem from a simple classical Hamiltonian. But is the usual spin-ice Hamiltonian a correct description of the experimental systems? Here we address this issue by measuring magnetic susceptibility in the two most studied spin-ice compounds, Dy2Ti2O7 and Ho2Ti2O7, using a vector magnet. Using these results, and guided by a theoretical analysis of possible distortions to the pyrochlore lattice, we construct an effective Hamiltonian and explore it using Monte Carlo simulations. We show how this Hamiltonian reproduces the experimental results, including the formation of a phase of intermediate polarization, and gives important information about the possible ground state of real spin-ice systems. Our work suggests an unusual situation in which distortions might contribute to the preservation rather than relief of the effects of frustration.

S pin-ice systems owe their name to an analogy of their ground-state construction rules to those in Pauling's model for proton disorder in water ice 1,2 . In their simplest variant, they can be described in terms of centre pointing classical Ising spins at the binding points of a network of corner-shared tetrahedra forming a pyrochlore lattice (see Fig. 1a). In real materials, such as Dy 2 Ti 2 O 7 (DTO) and Ho 2 Ti 2 O 7 (HTO), this is an accurate description for the low temperature behaviour. It was recognized early that, given the large magnetic moments of the rare-earth ions (of the order of 10 m B ) in addition to a ferromagnetic nearest neighbour coupling J (of the order of 1 K), a dipolar term was necessary in the Hamiltonian to describe the experimental results [3][4][5] . The minimal model that became the norm in the description of spin-ice materials is the standard dipolar spin-ice model (s-DSM): where r 1 is the nearest-neighbour distance, r ij is the distance between spins i and j, S i is a classical spin of unit length, and the dipolar constant, D, is of the same order than J. Spin-ice systems belong to the class of geometrically frustrated magnets. Frustration is possible in other ways, but the class of geometrically frustrated materials has the advantage of better experimental control and theoretical description than the case where disorder is needed. An inherent weakness of geometrical frustration is that it might be relieved by distortions to the lattice [6][7][8] . The spin-ice materials DTO and HTO are fairly robust to distortions; experimentally, the properties of the frustrated state and its low-lying excitations are seen to dominate the intermediate temperature regime 9 . However, as we shall see, the effects of small distortions in this system might account for some specific experimental features and determine the presence or absence of an ordered state at low temperatures. The different limits of the classical s-DSM already contain the most singular and attractive features of spin-ice materials, notably residual entropy, emergent gauge structure, fractionalized monopolar excitations [10][11][12] , and two-dimensional and threedimensional Kasteleyn transitions [13][14][15] . The s-DSM also predicts an ordered state at very low temperatures 16 . As the field developed, it was soon recognized that additional interaction  terms are needed to properly account for the experiments and an effort was made to determine an empirical Hamiltonian from the experimental data 17,18 . In particular, the most complete work of this type is that of Yavors'kii et al. 19 . In their approach, Monte Carlo simulations of a model with dipolar interactions, and isotropic first (J 1 ), second (J 2 ) and third (J 3 ) nearest-neighbour exchange interactions are compared with DTO experimental susceptibility, specific heat and magnetization data for zero magnetic field H 20-22 , with H applied in the 112 direction 23,24 , and with the empirical dependence of its polarization transition with H//111 25 . By means of this comparison they put bounds to the values of J 1 , J 2 and J 3 and arrive at an empirical Hamiltonian, the generalized dipolar spin-ice model (g-DSM). The g-DSM gives a very good phenomenological description of several of the spin-ice features and reproduces the low temperature neutron scattering pattern of DTO noticeably better than the s-DSM 19 . Despite its successes, the g-DSM model does not fully reproduce the low temperature behaviour of the spin-ice materials. In this work we address two instances of these, which were unknown at the time the model was developed: the doubling of the polarization transition with fields in the neighbourhood of 111 (refs 26,27) and the recent report of an upturn in the specific heat below 0.6 K 28 . Based on our experimental results and those in the literature, we construct a new spin-ice model and show that in addition to reproducing the experimental features, it predicts an intermediate magnetization state for fields near [111] and the existence of different types of possible order at low temperatures when no field is applied.

Results
Magnetic susceptibility. We start by discussing the magnetic susceptibility at low temperatures when an external constant field is applied in the neighbourhood of the crystallographic [111] direction. For our experiments, single crystalline samples of DTO and HTO were used. The a.c. susceptibility was measured in a dilution refrigerator with an external field applied using a tripleaxis vector (see Methods section).
The change in a.c. susceptibility at the polarization transition is plotted in Fig. 1 for the HTO single crystal (analogous behaviour is seen in the DTO sample) as a function of field at a fixed temperature of T ¼ 0.18 K for different angles from the [111] direction (see Fig. 1b). Figure 1c shows field rotations from [111] towards [110], while Fig. 1d shows rotations towards [112]. The usual spin-ice models predict a single transition between the partially frustrated kagome-ice state at low fields and the polarized state at high fields. In the experiments, this is true when the field is perfectly aligned with [111], or when it is rotated towards [110], but the transition splits into two peaks, with an intermediate polarization region, when the field is rotated towards [112]. This behaviour is seen more clearly in Fig. 2, where the change in susceptibility is shown for both materials as an interpolated contour plot. In this figure, it is also quite noticeable that the angular dependence of the lower transition field is quite small when the field is rotated towards [112]. Similar behaviour had been previously observed in magnetization experiments in DTO 26 .
The failure of the usual spin-ice Hamiltonians to reproduce an experimental feature common to both spin-ice materials suggests that there is an additional physical mechanism that is being neglected in the theoretical treatment. In the following, we argue that a missing ingredient can be the effect of distortions in the magnetic couplings of the pyrochlore lattice.
The role of distortions. We analyse how the s-DSM Hamiltonian (equation (1)) is modified in the presence of distortions. We consider the simplest case: classical Einstein phonon modes in the pyrochlore lattice, parameterized as u i Er i À r i 0 , with a linear dependence of the exchange couplings on distances: where J 0 is the undistorted exchange constant and a is the spin-phonon coupling constant, and a simple quadratic term for the energy cost of distortions. Expanding the variation of H up to linear order in u i , one obtains a magnetoelastic Hamiltonian that can be written compactly as: where K is the elastic constant for classical Einstein phonons, N(i) stands for the neighbours of site i and F ij encodes the quadratic spin interactions between a site i and a neighbour j. An effective Hamiltonian can be obtained from H me by integrating out the phononic d.f. 29  effect of distortions to the lattice is translated into a change in the nearest-neighbour interaction and the emergence of further neighbour interactions. These interactions have their origin in the local changes in the magnetic interaction due to changes in the relative positions of the lattice sites, and are thus strongly dependent on the lattice connectivity between neighbours. Further neighbours of equal distance but connected to the original site by a different number of lattice sites develop different interaction terms. In the case of the pyrochlore lattice (see Fig. 1a), this mechanism provides two kinds of third nearest neighbour interactions, J 3 and J 3 0 .
The model. While it cannot be expected that this simple analysis based on classical phonons will give accurate values for the effective exchange constants of the real system, which consists of a far more complicated lattice once all the non-magnetic atoms are taken into account, it provides the key observation that, regardless the specific material, further neighbour interactions should be considered and distinction made between differently connected neighbours of the same order, even if these terms or differences were negligible from the exchange integrals. With this information, we constructed a model to third nearest neighbours that incorporates these distinctions, the distorted dipolar spin-ice model (d-DSM), and explored it using Monte Carlo simulations (see Methods section and Supplementary Note 1). As additional constraints, we keep J 2 and the mean value (J 3 þ J 3 0 )/2, within the first set of intervals determined by Yavorsk'ii et al. for the g-DSM 19 to retain compatibility with the experimental results this model was fitted to (see also Supplementary Note 2 and Supplementary Fig. 1). Figure 3 shows the magnetic susceptibility for the d-DSM, with . 3a), and for the g-DSM model (Fig. 3b), both as interpolated contour plots. Unlike the s-DSM or g-DSM, and in keeping with the experiments, the susceptibility calculated in the d-DSM has a flatter angular dependence of the critical field for rotations towards [112], showing a broad peak that eventually resolves into two peaks, and remains unique when the field is rotated towards [110], all of which is in correspondence with the behaviour shown in Fig. 2 for the experimental results in the neighbourhood of [111]. Figure 3b shows the best possible fit to the data using the g-DSM: by forcing the parameters towards the higher end of the intervals given in ref. 19, it is possible to eventually achieve a small doubling of the transition when rotating towards [112], but it completely fails to reproduce the angular dependence of the critical field. The addition in the d-DSM of different types of J 3 is therefore necessary to be able to reproduce the susceptibility experiments close to [111], even at a qualitative level. Further terms and analysis are required to reproduce in detail the full angular dependence.
Intermediate phase. A known consequence of the coupling between spin d.f. and distortions is the stabilization of plateaux in the magnetization (see for example, refs 30,31). It therefore comes as little surprise that the doubling of the polarization transition corresponds, in the limit of low temperature and perfectly homogeneous tensions, to the presence of a plateau between the kagome-ice and the fully polarized state (see Fig. 4a). This plateau marks an intermediate ordered state where the ice rule is broken in only half of the tetrahedra. To visualize these states, it is helpful to mark by a red or blue dot the defect or excess of inward pointing spins. Defects of different colour interact approximately as charges of opposite sign 11 . Full polarization corresponds in this picture to a complete checkerboard pattern of blue and red dots. The intermediate states of Fig. 4 are easily pictured as alternate arrangements of stripes of the checkerboard pattern separated by empty stripes (see Fig. 4).
The ground state. One of the remarkable predictions of all dipolar spin models is that, contrary to initial expectations, spinice materials may have no extensive residual entropy even in the absence of an applied magnetic field. As shown by Melko et al. 16 , the introduction of dipolar interactions leads to a sharp transition at B0.18 K into an ordered ground state. At present there is no direct experimental determination of the ground state of spin-ice systems, since below 0.6 K the characteristic relaxation time of the system seems to grow faster than exponentially 32,33 and most measurements probe properties of the system out of equilibrium 34 . A recent experiment 28 , in which the zero field low temperature specific heat was measured over much longer timescales, reports the possible observation of a peak in the specific heat that might correspond to the onset of order in DTO. However, the temperature at which this onset is seen (B0.6 K) is not compatible with the prediction of Melko et al. for the s-DSM. The situation is very similar for the g-DSM: the additional interaction terms of this model leave the ordered state unaltered, and the specific heat (by construction), shows no upturn down to 0.2 K.
The model presented in this work (the d-DSM) is a different scenario altogether. In the simple homogeneous case, and within the constraints on the parameters imposed by the experiments (crucially our susceptibility near [111]) two different ordered ground states are possible depending on the relation between the exchange constants. Figure 5a,b shows cuts along [100] of the two ordered states obtained at low temperatures for different values of the exchange constants. For small values of J 2 , the relevant parameter that switches between the two types of order is the difference, The scale in the middle (Fig. 5c) indicates, within the region compatible with the experiments (shaded in blue), the values of DJ 3 (for fixed average J 3 , and J 2 ¼ 0) where each type of order is seen. As expected, for low values of DJ 3 , the ordered state is identical to the one originally found in the s-DSM. A simple description of this state (state I) is to consider it an antiferromagnetic arrangement of ferromagnetic chains of spins (indicated with the same colour in the figure) at axes rotated p/4 with respect to the original lattice. The second state (state II) can be thought as a simple extension of state I where the unit blocks are now pairs of ferromagnetic spin chains (also indicated with the same colour in the right hand side figure) ordered antiferromagnetically. In both cases, the non-shaded region in the figure corresponds to the unit cell.
These are the two states of lower energy that have been identified depending on the degree of asymmetry in the two J 3 . It is straightforward to conjecture that there is a wealth of possible intermediate states with energies differing by a few mK (in 1/k), constructed from all the possible periodic intercalations of different numbers of alternating single and double chains in each direction, or even of the alternation of bundles of a higher number of chains. This by itself has already the consequence that it will be extremely difficult for any experiment to probe equilibrium properties, but the situation in a real sample is even more complicated: the most likely experimental state is that the sample will have a distribution of tensions and distortions within its volume. This will lead to local changes in DJ 3 , coming from the distortion-induced contribution to the two J 3 , which in turn will result into local energetic rearrangements of the possible states. No long-range ordered state will be favourable in the whole sample. The effect of tensions and distortions could also be intimately related with the experimental observation of freezing or glassy behaviour below temperatures of about 0.6 K 35-38 .
Specific heat. Figure 5d shows the zero-field specific heat divided by temperature, calculated for the d-DSM with the same values used for Fig. 3, compared with the experimental data for DTO from refs 12,28 (a finite size analysis of the simulated specific heat is shown in Supplementary Note 3 and Supplementary Fig. 2). Both experimental sets of data coincide at high temperatures, down to about 0.6 K below which the characteristic time of the system becomes extremely long and the results from ref. 28, taken over longer times, are believed to be closer to equilibrium. The upturn seen in this data has been interpreted as a possible sign of onset of long-range order. The specific heat calculated from the d-DSM matches the high temperature behaviour of the experiments and shows a sharp peak at a temperature close to 0.2 K, which corresponds to the ordering transition into state I or II (depending on the value of DJ 3 ). It would be tempting to associate this peak with the upturn of ref. 28, and match the two sets of data by making small adjustments to the exchange constants of the d-DSM. However, as discussed before, it might not be realistic to expect a homogeneous long-range ordered state in a real sample, but rather a wide range of local configurations and a multitude of energetically close excitations. In this case, features matching those seen in the specific heat data of ref. 28 have their origin in a Schottky-type peak consequence of this proliferation of states of very similar energies in this range of temperatures. A more definitive answer to this question would come from neutron scattering experiments performed under the same relaxation conditions as the specific-heat experiment.

Discussion
In summary, in this work we present experiments performed in the two paradigmatic examples of spin-ice materials, DTO and HTO, where we have measured the angular dependence of susceptibility using precise alignment with the crystallographic axes, concentrating in particular in the polarization transition with field in the neighbourhood of [111]. We show that a generic feature of these two materials, absent from any theoretical spinice model in the literature, is a marked anisotropy in w, which shows a doubled transition when the field is rotated towards [112]. This doubled transition had been previously observed in specific heat and magnetization experiments in DTO 23,26 , but not in HTO. We considered the possibility of distortions, a missing ingredient from the usual spin-ice models, and by performing the simplest possible analysis on their effect on the s-DSM we found that they provide a justification for the addition of interactions between second and two kinds of third nearest neighbours. This is a generic feature of distorted spin-ice materials, and applies even if the exchange terms for the specific material were negligible. There is a crucial difference in the current treatment of thermal vibrations to what already exists in the literature, which is a consequence of the presence of long-range interactions in this class of materials. Thermal vibrations usually just renormalize constants, but in the present case dipolar long-range interactions lead to the appearance of couplings that would not exist in the absence of vibrations. Thus, our treatment of the simple d-DSM leads to a more complex Hamiltonian, that resembles that of for example, ref. 19, but is not obtained in a purely empirical manner.
Provided with these ingredients we constructed a new model (the d-DSM) and explored its characteristics using Monte Carlo simulations. The new model reproduces the anisotropy found in the experiments with fields in the neighbourhood of [111], and, despite the fact that it was not built based on neutron scattering data or Cv, the S(q) calculated for the d-DSM reproduces the experimental results (see Supplementary Note 4 and Supplementary Fig. 3) and we have seen that the same is true for Cv (Fig. 5). More importantly, it predicts an intermediate magnetization state, thus explaining the hitherto unresolved issue of a doubled polarization transition with field close to the [111] direction, and also predicts the existence of different types of possible order at low temperatures when no field is applied. These results from the d-DSM are valid regardless the ultimate origin of the different interaction terms. The possibility that these further order terms could be a consequence of distortion would mean that the small tensions present in real samples could result in the absence of long-range order at low temperatures. This is probably the most surprising conclusion of this work, that distortions, which usually result in a relief of frustration, might be the reason why the system would remain disordered at very low temperatures. Very recently, Henelius et al. 39 have presented an empirical spin-ice Hamiltonian constructed to describe three main sets of experiments in DTO (the structure factor S(q), the specific heat Cv and the magnetization for H near [112]). The model Hamiltonians presented in both studies have several points of contact. In particular, both models predict the same possible ordered ground states for H ¼ 0 (those described in Fig. 5) and, crucially, the likelihood of a disordered ground state as a consequence of internal stress (in our case given by external strain, in theirs by impurities) and quasi-degeneracy. The similarity of these findings is remarkable taking into account that the core experiments used in their work and in ours are very different.

Methods
Experiments. For our experiments, single crystalline samples of DTO and HTO where grown using an optical floating zone in St. Andrews and Oxford, and they were oriented and cut into prisms of approximate dimensions (0.7 Â 0.7 Â 3) mm 3 , with the long axis oriented along [111]. The a.c. susceptibility was measured in a dilution refrigerator using two sets of detection coils, positioned so as to avoid cross-talk effects, and a drive coil. Each detection coil set comprises of two coils, connected antiparallel to each other. The drive coil is concentrically wound around the pick-up coils. The samples were placed in each of the upper pick-up coils and thermally grounded to the mixing chamber through silver and copper wires. A low-frequency and very low-amplitude excitation field of about 3 Â 10 À 5 T r.m.s. was generated by the ac current in the drive coil, and the response detected by a lock-in amplifier. For accurate control of the external magnetic field, we used a bespoke triple-axis vector magnet capable of producing 9/1/1 T along z, x and y. Careful centring using the vector field capability allowed y to be determined to an accuracy of 0.1°relative to the crystallographic axis [111].
Theoretical analysis of the effect of distortions. A simpler effective Hamiltonian can be obtained from H me (equation (2)) by integrating out the phononic d.f. 26 . In the spin-ice model S i ¼ s i e i , with s ¼ ± 1, which leads to important simplifications: since s i 2 ¼ 1, the usual quartic terms become quadratic. Still, due to the long-range nature of dipolar couplings one encounters terms F ij of any range, and the analysis of the effective magnetic model is rather complex. To compare the predictions of our model with those in the literature we have considered an expansion of the corrections to the full effective magnetic Hamiltonian up to third neighbours. The effective magnetic model obtained reads where oi,j4 r indicates r-th range neighbours. The couplings up to third neighbours are where d is the contribution arising from the phonons, d ¼ J=4ka 2 ða a=3 À 5D=JÞ 2 . Further range terms are the unmodified dipolar interactions. The first observation is that the crystal environment distinguishes two different third neighbour couplings in H eff : J 3 eff and J 3 0eff distinguish between third nearest neighbours connected through the lattice by a minimum of one or two other lattice sites respectively. One of them is sensitive to distortions, while the other is not. A second observation is that the relative strengths of all couplings in H eff depend only on two microscopic parameters: D/J and a/J. Corrections of longer-range magnetic couplings can be done systematically, depending on the same microscopic parameters. We have checked that fourth-order corrections do not change the main results of this work.
It is not realistic to expect that a simple classical model such as the one presented here will give an accurate description of the real system at low temperatures. The presence of soft modes (speculated for spin-ice materials in ref. 40 where the oi,j4 3 and oi,j4 3' distinguish between third nearest neighbours connected through the lattice by a minimum of one or two other lattice sites respectively (see Fig. 1a). We performed Monte Carlo simulations, using Ewald summations to take into account long-range dipolar interactions 42,43 in a conventional cubic cell for the pyrochlore lattice. We simulated systems with L Â L Â L cells with periodic boundary conditions. Basic parameters were taken from ref. 19, including cell size, lattice parameter and magnitude of magnetic moments. Thermodynamic data as a function of field were collected using a variation of single spin-flip Metropolis algorithm. For Fig. 3 we chose a temperature of 0.4 K to avoid problems associated with hysteresis. For the specific heat data below 0.6 K we mixed the usual Metropolis with a dynamics that conserves defects for a certain number of Monte Carlo steps before letting the system relax (the algorithm is a variation of the one used in ref. 43 and is described in Supplementary Note 1). In the case of magnetization measurements at 0.1 K (Fig. 4), data at each field were collected after a slow annealing, averaging over five independent runs.
The average value of the third neighbour interactions was kept within the range specified in ref. 19 to retain compatibility with previous experiments in DTO. As an example, with the values of J used in Figs 3-5, we studied the transition to ferromagnetic ordering with field nearly along the [112] direction 19,21 . We obtained a transition temperature of E0.35 K, which could be justified considering angular deviations of o0.5°off perfect alignment.
Regarding the double peak near [111], we contrasted our simulations with the static magnetization data obtained in ref. 26, which shows narrower features than the a.c. susceptibility. The range fixed in Fig. 5 was determined looking for special key details: (i) a magnetization increase of only B0.025m when the field is tilted from 0°to À 5°(towards [112]) at a fixed field of 0.9 T, (ii) an increment for the first critical field of o0.06 T when tilting the field from 0°to À 5°and (iii) a difference between the two critical fields at À 5°of E0.1 T.
Data availability. All relevant data are available from the authors.