Monopole matter from magnetoelastic coupling in the Ising pyrochlore

Ising models on a pyrochlore oxide lattice are usually associated with spin ice materials and"magnetic monopoles". Ever more often effects connecting magnetic and elastic degrees of freedom are reported on these and other related frustrated materials. Here we extend a spin-ice Hamiltonian to include coupling between spins and the O$^{-2}$ ions mediating superexchange; we call it the Magnetoelastic Spin Ice model (MeSI). There has been a long search for a model in which monopoles would spontaneously become the building blocks of new ground-states: the MeSI Hamiltonian is such a model. In spite of its simplicity and classical approach, it describes (both spin and oxygen lattice) the double-layered monopole crystal observed in Tb$_2$Ti$_2$O$_7$. Remarkably, the dipolar electric moment of single monopoles emerges as a probe for magnetism. As an example we show that, in principle, pinch points related to Coulomb phases could be detected in association with the O$^{-2}$-ion displacements.


I. INTRODUCTION
It is remarkable that even one of the simplest interacting models in condensed matter physics (the Ising model) can lead to some of the phenomena in geometrically frustrated magnetism that has kept busy part of our community for the past decades. At the core of the new emergent physics -massively degenerate ground states with critical-like spin correlations, exotic excitations, artificial electrostatics, and very peculiar dynamics [1-7]-is the strategic choice of the spin lattice structure so that pairwise interactions compete rather than collaborate.
The pyrochlore structure (Fig. 1a)) is a prominent example among these "frustrated" lattices, with spin ice canonical materials Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 as some of its most notable examples [5][6][7][8][9][10]. Their effective residual magnetic entropy is similar to that of water ice, and the source of their collective name. The configurations of Ising spins in the lowest energy states of these materials [11][12][13] can be described by a lattice gauge field that fluctuates like an electric field in vacuum [4,[14][15][16][17]. The combination of this "Coulomb Phase" with nonnegligible dipolar interactions leads in turn to the emergence of local magnetic excitations: the "monopoles". They sit in the centres of the tetrahedra that make the pyrochlore lattice, and interact through Coulomb forces like electrical charges [12,18]. As illustrated in Fig. 1a) there are different types of these magnetic charge-like quasiparticles (eight "single" monopoles, two "double" ones). Monopoles are responsible for the very peculiar dynamics of spin ices at low temperatures [19][20][21][22]. Also, under different conditions, they can act as building blocks * These authors contributed equally.
† Corresponding author e-mail: borzi@fisica.unlp.edu.ar for different "monopole phases" [8,12,[23][24][25] that have been studied theoretically or experimentally. In general, dense "monopole matter" was forced to appear by resourcing to somewhat artificial conditions [23,[26][27][28][29], freezing spin fluctuations [12,24], imposing out of equilibrium situations [30], or breaking some symmetry of the system [28,[31][32][33][34][35]. It can be proved to be impossible to obtain the most general monopole liquid solely from pairwise interactions [29], leaving unanswered a fundamental question that we pursuit to respond here: how can monopole matter be thermodynamically stable in real materials without explicitly breaking any symmetry? Central to this question and to this work is the interplay between magnetic and elastic degrees of freedom. Since it is the precise geometry of the lattice the one that balances out the pairwise spin interactions, geometrically frustrated systems can be quite susceptible to spontaneous deformation [36][37][38][39][40][41][42]. Regarding Ising pyrochlores that remain disordered at the lowest temperatures, this coupling is responsible for structural fluctuations [43], giant magnetostriction [44,45], and composite magnetoelastic excitations in Tb 2 Ti 2 O 7 [46]. It seems to be much smaller in the canonical spin ices [47,48], but may explain subtle effects shaping the zero magnetic field (h = 0), and h||[111] phase diagrams of Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 [49], dynamics [50], and the observed magnetic avalanches [21,51,52]. Khomskii [53] was the first to notice that, in spin ice, spin configurations related to single monopoles are necessarily accompanied by local distortions that result in an electric dipole. These dipoles can interact with an external electric field [53] or among themselves [54], changing the energy balance.
Inverting Khomskii's line of reasoning, we demonstrate in this work that magnetoelasticity can be the keystone for monopole stabilization in pyrochlore oxides. We modify the nearest neighbours spin ice Hamiltonian in the simplest possible way to include a coupling to the lattice of O −2 -ions sitting near the centre of tetrahedra (see Fig. 1). In the regime of strong coupling, lattice distortions turn the eight types of single monopoles into stable, atomic-like constituents of novel ground states. In our model, distortions are not just dummy variables but dynamic degrees of freedom. We can contrast their behaviour with that of real materials, employ them as probes to investigate the underlying magnetism, or -in a multiferroic fashion-to control the material properties using electric fields.
In spite of its simplicity, and building on the previous works of Jaubert and Moessner [54] and Sazonov and collaborators [25], our Magnetoelastic Spin Ice (MeSI) model allows to understand in a new manner the formation of a double-layered monopole crystal in Tb 2 Ti 2 O 7 with field applied along [110], and to contrast the O −2 -distortions with those suggested previously [25]. The model's output is compatible with the power-law spin correlations observed in Tb 2 Ti 2 O 7 at zero field [55], and gives some clues on the half moons measured in neutron diffraction patterns at finite energy [46,56].
Although we concentrate on the strong coupling limit, we expect the MeSI model to be a convenient tool to study other systems, in particular spin ices. Incorporating the lattice degrees of freedom may open the way to the survey or design of the electrical properties on Ising pyrochlores, or teach us how to probe other properties through them (as it has been done in some pioneering works in spin-ice [57][58][59][60]).
The article is organised as follows. We begin the Results Section by introducing and justifying the extended magnetoelastic model (Subsection II A). We then show how the MeSI model stabilizes a Monopole Liquid (Subsection II B) . This massively degenerate perfect paramagnet will serve as the basis from which the other cases of study will follow through small perturbations. Including attraction between monopoles of equal charge will lead to a phase comparable to the "jellyfish" or "spin slush" [30,61], with half moons in the neutron structure factor. Correspondingly, Coulomb-like attraction gives rise to a Zincblende Monopole Crystal with magnetic moment fragmentation [26,28] (Subsection II C) . We will see that the deformed O lattice fluctuates with the fragmented magnetic moments; the pinch points in its structure factor would be, in principle, detectable using diffuse x-ray diffraction. The explicit breaking of a symmetry by a magnetic field is addressed in Subsection II D (the double-layered Monopole Crystal [25,54]). The Discussion (Section III) shows some striking results concerning the structure factors of different phases, evaluates the possibility of observing these effects, and suggests possible new avenues of research.

A. A model for magnetoelasticity in Ising pyrochlores
We will study a pyrochlore oxide lattice with Ising spins of the type R 2 M 2 O 7 . Spins will generally be associated with rare earth ions (R), while M is typically a transition metal (Ti, Sn, Zr) [62][63][64] but could also be Ge or Si [65]. The spins sit in the corners i = 1...4 of up-tetrahedra (coloured purple in Fig. 1a)). They point along the 111 directions, either towards (with pseudospin variable S i = 1) or against (S i = −1) the centre of the tetrahedron they belong to. In order to better describe the variety of states of matter we are going to study, it will be useful to employ the language of magnetic excitations or "monopoles". We will use these two terms to refer to the topological charge even in the absence of long range dipolar interactions [12]. One can group the different spin configurations of a single tetrahedron into sets, using the net entrant spin flux as a label that defines the magnetic charge [12] of that tetrahedron. The same definition is valid for up or down tetrahedra. A crucial observation is that fixing the magnetic charge in a tetrahedron does not necessarily define the spins variables in a unique way. There are six different "neutral" or "spin ice" configurations, with two spins pointing in and two pointing out (empty tetrahedra in Fig. 1a)). There are four positive (negative) single monopoles of charge Q (−Q), with three spins pointing in and one out (three out and one in); they are drawn as small green (red) spheres in Fig. 1a). Finally, for double monopoles each charge identifies a single configuration: 2Q (−2Q) when all the spins point in (out) of the tetrahedron; a negative double monopole is represented by a big red sphere in Fig. 1a).
Each tetrahedron in the pyrochlore lattice can be embedded in a cube. The six links between nearest neighbour spins lie diagonally along the six faces of the cube and can be labelled using the perpendicular Cartesian axes (e.g. +z and −z for the links between spins S 1 -S 3 and S 2 -S 4 , shown in green and red respectively in Fig. 1b)). Following other studies [54,[62][63][64] we assume that the superexchange between R-ions takes place through the oxygen ion O −2 , sitting at the centre of the tetrahedra (see Fig. 1b)). In order to simplify our model for magnetoelastic coupling we will only consider the independent displacement of these non-magnetic ions, keeping all the rest at fixed positions. The restoring force for the oxygen points towards the centre of the tetrahedron and is taken to be isotropic and proportional to the oxygen's relative displacement δr (see Fig. 1). With these considerations, and taking into account only nearest neighbor magnetic interactions, our model Hamiltonian can be written as Here δu ≡ δr/r nn , with r nn the nearest neighbor distance, K is the elastic constant for the oxygen ions. The sum runs over all (up and down) tetrahedra. J ij (δu) is the displacement-dependent nearest neighbors superexchange energy associated to each pair. It can also be labeled using the link name J ±m (δu), with m = x, y, z (for example, J +z ≡ J 13 for the up tetrahedron in Fig. 1b); for more details on the notation see Supplementary Information I). For small deviations δu, the superexchange constants can be expanded around the undistorted value [42], J 0 , which corresponds to the configuration where the O occupies the central position and is thus identical for all directions. We will assume that the main effect of the O displacement over the exchange constants comes through the change in the bond angle of the R-O-R [53,54]. The net result of the angular distortion on J is to make it more antiferromagnetic or ferromagnetic according to the Goodenough-Kanamori-Anderson rules. As shown in Supplementary  is the coupling constant of the global system that correlates the lattice and magnetic degrees of freedom.
Within this approximation it is possible to recast the Hamiltonian (see Supplementary Information I) into a compact vectorial form, where the dependence on the O −2 lattice distortion is explicit. We call this extension of the simplest spin ice model the Magnetoelastic Spin Ice model (MeSI); for zero magnetic field it is given by: Here we have defined J ml ≡ 3α 2 /K and δũ ≡αδu, both measured in Kelvin; the sum runs along the diamond lattice, J 0 = (J 0 , J 0 , J 0 ) and the vectors [S,S] ± have components The first term of Eq. 2 is the elastic energy, and it is easy to see that the last one is the usual nearestneighbour Hamiltonian with isotropic exchange constants. If different types of nearest neighbour bonds were to be considered (as we will do when considering the effect of magnetostriction in Section II D) the latter would be replaced by a sum involving the exchange constant matrix J ij 0 : i =j J ij 0 S i S j . The middle term in the MeSI model is central to this work, as it contains the (linearized) magnetoelastic coupling. While the coupling constant is somewhat hidden inside δũ =αδu, we will soon show that the new energy scale J ml (proportional to the square of the coupling constantα) is a convenient measure of the relative stability of single monopoles, the atomic-like building blocks of the new exotic phases we will study in the following sections. Also, and equally important, it indicates how strongly magnetism will be reflected in structural properties and measurements. Models of interacting entities, even simple ones, are seldomly exactly solvable. It is then a surprise that the three-dimensional MeSI model turns out to be analytically solvable for J 0 = 0. Completing squares in δũ, the Hamiltonian can be decomposed into an "elastic" and a "magnetic" term. The first one is where the components δO m = δũ m − J ml 3 [S,S] m − can be interpreted as the relative displacement of the oxygen with respect to its equilibrium position along the different axes. Due to the magnetoelastic coupling this position depends now on the specific spin configuration in each tetrahedron. This term is quadratic and can be easily integrated out. If we include a Zeeman term, proportional to the magnetic field h (measured in Kelvin), the effective magnetic term under a magnetic field then becomes The last two terms are the nearest neighbour spin-ice Hamiltonian under an applied magnetic field (with uniform exchange constant J 0 ) . For a strong magnetoelastic coupling, the first term (with a four-spin product) stabilizes a Monopole Liquid at low temperatures [29]. It is easy to check that the range of stability is given by J 0 < J ml for positive J 0 (which otherwise corresponds to a spin ice phase), and J ml > −3J 0 for negative J 0 (usually leading to a double monopole crystal).
The Monopole Liquid for J 0 = 0 has been shown to be a perfect paramagnet, with no spin correlations at any temperature [29]. Its ground state holds a massive residual entropy, and is equally populated by the 8 possible monopole configurations. The four-spin model (i.e., H eff for h = 0 and J 0 = 0) was solved exactly by Barry and Wu ten years before the discovery of Spin Ice [66]. In recent years it had been suggested the possibility that lattice distortions could stabilize dense monopole phases [25,28,29]; the MeSI model crystallises this idea in a clean and straightforward fashion, with the added benefit of an analytical solution. Fig. 2 shows results of our Monte Carlo simulations (symbols) for the full MeSI Hamiltonian for J 0 = 0 (Eq. 2). They are compared with the exact results obtained by Barry and Wu [66], displayed as full lines. Note that, unlike the case of Ref. 29, the model here involves both the spins and the (coupled) elastic degrees of freedom. We define the density of monopoles, ρ, as the average number of single monopoles per tetrahedron (note that it does not count double charges). Fig. 2 shows that ρ(T ) saturates at ρ = 1 monopole per tetrahedron for T /J ml << 1, as expected for a dense phase of single charges; on the other hand, the inverse magnetic susceptibility χ −1 is that of a paramagnet, with no evidence of increasing magnetic correlations with decreasing T (Fig. 2a)). In both cases there is very good agreement between the simulations of the full model and the analytical results for the effective model [66]. On the other hand, the specific heat per unit spin C V and mean square deviation δu 2 (Fig. 2b)) make apparent that we are in fact dealing with a composite magneto-elastic system. The solution of Barry and Wu for C V /k B needs an extra constant term of 3/4 to take into account the elastic energy of the N/2 oxygen ions, as expected from the equipartition theorem. Although according to this same theorem one would naively expect a straight line for δu 2 The cubic lattice size we simulated had N =8192 spins and N/2 moving oxygens; J ml = 3α 2 /K = 1K, to fix the energy scale. a) As noted in Ref. [29], the spin part behaves as a perfect paramagnet; there are no spin correlations at any temperature, with a perfect Curie-law for the inverse susceptibility χ −1 . The density of monopoles ρ saturates at 1 single monopole per tetrahedron for T << J ml (i.e., forms a "Monopole Liquid"). b) The specific heat per unit spin CV and the mean quadratic deviation δu 2 evidence that we are really dealing with a composed system. Due to the elastic contribution H el from the oxygen ions, we needed to add a constant term 3/4 to the analytical solution to match the simulations. δu 2 (T ) is not just a straight line, owing to the dependence of the O equilibrium position with the spin configuration in δO.
vs T , there is a kink for δu 2 in Fig. 2b) for T below the maximum in C V . It is a sign of coupling between the degrees of freedom: as we mentioned, the oxygen equilibrium position depends on the local spin configurations (Eq. 4).
The interplay between the nearest-neighbours spin ice term proportional to J 0 and the four-spins term favouring a monopole liquid has been explored in Ref. 29, where the four spin term in Eq. 5 was proposed as a model Hamiltonian. In addition to the Zeeman term included in Eq. 5 (studied in the next sections) it is interesting to consider interactions between monopoles (as those that would arise by including magnetic dipolar interactions between spins [12]). A simple way to introduce near-est neighbour repulsive or even attractive forces between like-monopoles consists in including second and third nearest neighbours spin interactions with a carefully chosen ratio [30,67]. In order to preserve generality, we will express the interaction directly in terms of the monopolar charge on a tetrahedron, Q : where Q = 0, ±1, ±2. Here the sum runs over nearest neighbour tetrahedra, and γ measures the strength of like-charge repulsion (γ > 0) or attraction (γ < 0). We have referred to the liquid phase with a single monopole per tetrahedron (ρ = 1) and no spin correlations as "the" Monopole Liquid, ML. However, other monopole liquids can be obtained by applying an external magnetic field [35], changing the monopole composition, or the spin or charge correlations [29]. Some of these phases show particular patterns in the structure factor that have led to different monikers. "Half moons" or "split rings" have been observed in the structure factor at the "jellyfish point" [30] or the "spin slush" phase [61], with single monopole density ρ ≈ 0.35 and attraction between like-monopoles. We have calculated the neutron structure factor I Spin (k) (see Supplementary Information III for details) for the ML in the presence of nearest neighbour attraction between like charges as per Eq. 6 with γ < 0 (ρ = 1, T << |γ| < J ml ). The diffuse pattern we obtain can be understood as the merging of the different half moons observed in Refs. 30 and 61, as their features widen due to the higher density. While half moons are usually detected at finite energy [56,68], the ML with like-attraction is a new instance (together with Refs. 30, 61, 69, 70) of a ground state with this feature. Interestingly, the need for an attraction between like monopoles will arise again when studying Tb 2 Ti 2 O 7 (Section II D), a compound which also shows half moons in its neutron structure factor [46,56].
With their exotic excitations [12,18], topological phase transitions [7,16,[71][72][73][74], peculiar dynamics [3,[19][20][21][22]75], power law correlations leading to pinch points [4,14,16,17,76,77], the possibility to tune new ordered or disordered phases using magnetic fields [7,8,13,49,[78][79][80][81][82][83][84], spin ices have showed a wealth of interesting physics. In the same way, the opportunity to stabilize a completely different phase with massive residual entropy in an Ising pyrochlore opens the door to new and non-trivial forms of dense monopole matter. Part of these phases have been theoretically speculated on [12, 23, 26, 28-30, 32, 33, 35, 54] or inferred through experiments [25,31,49]. In what follows, we analyze how to obtain from the MeSI model some of these states, which, even within the realm of classical systems, do not exhaust all the possibilities opened by the inclusion of the magnetoelastic coupling. There exist previous experimental realisations and theoretical proposals for the single Monopole Crystal with the Zincblende structure (ZnMC, for short), stabilised at low temperatures by means of extrinsic [33-35, 78, 85] or internal fields [26,28,31,86], or dynamical constraints [23]. Within this context, it was first established that spins in a crystal of single monopoles at zero field could still fluctuate [23]. Brooks-Bartlett and collaborators noted that these partially ordered spins fragmented into two independent parts [26]. A static divergence-full part was related to the monople crystal, and the remaining (divergence-less) fragment, with neutron pinch points [26,28], characterised a Coulomb phase [4]. A number of Fragmented Coulomb Spin Liquids (FCSL) have been recently achieved experimentally in a pyrochlore lattice [31,86]. There, the Ir sublattice orders antiferromagnetically, acting as an effective field (with staggered values on up and down tetrahedra) over the spins in the other pyrochlore sublattice (Ho and Dy, respectively).
Returning to our work, the inclusion of an effective monopole attraction between + and − charges (γ > 0 in Eq. 6), implicit, for example, on dipolar spin interactions, will transform the fluid of single monopoles studied in Section II B into a ZnMC on decreasing temperature. As studied before [26], this phase would show magnetic moment fragmentation, with pinch points in the diffuse structure factor. However, in contrast with previous cases, there should now be spontaneous symmetry breaking between the two sites of the diamond lattice. The staggered charge density ρ S (defined as the modulus of the total magnetic charge due to single monopoles in up tetrahedra per sublattice site per unit charge) is the order parameter of the transition, which has a complex phase diagram [23,87]. Fig. 3 shows Monte Carlo simulations for the MeSI model with J 0 = 0 and opposite sign attraction (γ/J ml = 0.2) in Eq. 6. We observe a high density of monopoles in the whole temperature range, saturating at ρ = 1 at low T . The peak in the specific heat reflects the formation of a crystal (panel a)). Contrary to previous studies, [31,32], the abrupt increase in ρ S , together with the peak in its fluctuations shows that this time the symmetry between up and down tetrahedra is spontaneously broken at the transition ( Fig. 3 panel b)). By varying the value of J 0 > 0 the whole phase diagram ρ vs. T for charges in a lattice obtained for electric charges in a lattice [88], and then for conserved magnetic monopoles [23] can now be understood as emerging from a classical Hamiltonian with physical foundations. Including a negative J 0 , a double monopole crystal can also be stabilized [24,89]. Defining the total density of monopoles ρ T to include double monopoles (0 ≤ ρ T ≤ 2), a complex phase diagram would then be obtained comprising three different ground states: the vacuum of monopoles with ρ T = 0, the crystal of single monopoles for ρ T = 1 (both exponentially degenerate and with an associated gauge field, if no other interactions are added), and the zero-entropy crystal of double monopoles for ρ T = 2.
Even if we take it as a possible route to the relatively unexplored physics of "condensed monopole matter", we would not be making justice to the MeSI model if we do not consider in more detail the new, structural degrees of freedom. As we will see below, this allows us to make apparent some of the consequences of fragmentation from a different perspective. As sketched in the inset to Fig. 3b), when monopoles are stabilized at low temperatures the oxygen ions tend to be displaced along the 111 directions, towards one of the four triangular faces of the tetrahedron. This is the face that contains the three antiferromagnetic-like links (out-out, or in-in), painted green in Figs. 3b) and 4b). It is easy to check that the O −2 displacements δu i of a positive single monopole points antiparallel (parallel) to the total magnetic moment µ i of an up (down) tetrahedron. On reversing time the magnetic charge and dipole moment invert their direction, but the displacement δu i (and the electric dipole moment) remains fixed. For a given monopole charge, then, electric and magnetic moments flip in unison.
We can then argue that since the divergence-less part of the magnetic moment fluctuates like a gauge field, with neutron scattering pinch points in its structure factor [26,31], the same should be true for the dipolar electric moment sitting in each tetrahedron. If this were true, aside from the usual Bragg peaks associated with the pyrochlore crystal, pinch points related to correlated oxygen fluctuations around the centre of the tetrahedron could in principle be detected as diffuse scattering using simply an electron beam, or x-ray diffraction. The chances to observe the effect depends critically on the magnitude of the O −2 displacement. We will discuss this in Section III, where we also summarize the structure factor results for this and two other phases.
The fact that electronic dipolar magnetic moments could give birth to magnetic charges, and that these magnetic entities have associated electric dipolar moments has been mentioned as a further remarkable example of symmetry between electric and magnetic charges [53]. Another layer of complexity is thus added by noting that the correlated fluctuations of these dipolar electric and magnetic moments lead to twin gauge fields, that could be measured by probes coupling either to electric charge or to magnetic moments.  [25,90] (see Fig. 4a)). To justify this charge order, Jaubert and The order parameter for crystal formation (the staggered monopole density, ρS) and its fluctuations χS in the same temperature range. We have subtracted the contribution of the pure vibrational degrees of freedom from CV . The peak in χS reflects the spontaneous symmetry breaking. The inset to panel b) shows a magnetoelastic configuration of minimum energy for a positive single monopole in an up tetrahedron. The total magnetic moment (thin yellow arrow) points along the minority spin (blue arrow); it is one among a total of 8 different directions: 4 associated with positive and 4 with negative single monopoles, or 1 per spin configuration. The short black arrow corresponds to δr, the O −2 displacement of minimum energy for this magnetic configuration. There are only 4 different displacements δr, since they do not invert under time reversal, and they are always directed towards the face of the tetrahedron with 3 antiferromagnetic (i.e. in-out) links, coloured green here.
Moessner [54] explored a classical model. The mechanism involved the long range interactions between the electric dipoles associated with single magnetic monopoles [53], and in a much lesser degree magnetic dipolar interactions. They found a transition from the antiferromagnetic "allin/all-out" phase into the bi-layer when applying a [110] magnetic field.
Our MeSI model constitutes an alternative to this first intrinsic mechanism ever proposed to stabilize a single monopole phase [54]. While it can obviously not take into account all the complexity observed in Tb 2 Ti 2 O 7 [55,64,[91][92][93][94], it is an improvement over this previous proposal, since now both magnetic and elastic degrees of freedom are considered in the same footing. Furthermore, the model provides a unified explanation for the ground state observed at zero field (a Coulomb phase [55,95]), the double layered monopole crystal measured at moderate fields [25] (correcting the previously proposed O −2 lattice distortions), and suggests a connection with the presence of "half moons" in neutrons diffuse scattering at finite energy [46,56].
The application of a strong magnetic field along [110] does not fully order the Ising pyrochlore lattice. Spins on α-chains -running along [110], represented by purple arrows in Fig. 4a)-are completely polarized at low temperature. This results in effectively decoupled spins in β-chains (yellow arrows in the figure), with magnetic moments perpendicular to h. Only four possible spin configurations are then possible in each tetrahedron. Two of these are spin ice-like, with no average O displacement [79]; the other two are a positive and a negative single monopole, leading to the antiferromagnetic β−chains of spins, and β Q −chains of alternating charge shown in Fig. 4a) [33,34]. Unlike the figure (chosen to show a double monopole layer ordering) these β−chains are not coupled by the spin structure and -unless an explicit energetic coupling is included-would lead to an incoherent arrangement of β Q −chains.
Before tackling the question of charge coherence, the non-trivial question we should answer concerns magnetic charge stability. Why would Tb 2 Ti 2 O 7 change its ground state under a magnetic field h||[110] from a subset of the 2-in 2-out manifold to that of a dense phase of single monopoles [25] when the field is strong enough? Since the component of the magnetic moments along h is the same for the two chosen single monopoles and for the neutral sites, the response cannot rely on the Zeeman energy alone. It is interesting to note firstly that if we impose the alternate oxygen displacements along z-axis proposed by Sazonov et al (Fig. 6 in Ref. 25), the MeSI model naturally leads to a dense phase of non-coherent β Q −chains of magnetic monopoles. The only requisite is a displacement big enough to overcome the energy associated to the usual nearest neighbors term proportional to J 0 . Alternatively we will now inquire on the effect of h on the lattice structure, and then of that on the spin lattice through magnetoelastic couplingα.
We are not the first to notice the possible importance of the giant magnetostriction observed for h||[110] [45] for stabilizing the double-layered monopole phase in Tb 2 Ti 2 O 7 [54]. Here we will include its effect implicitly in the MeSI model through the exchange matrix J ij 0 in Eq. 2. Adding a Zeeman term for h||[110] the extended MeSI model can be written as: Rare-earth ions usually have associated a very strong spin-orbit coupling; through it, the torque acting on spins can affect the orbital angular momentum, and then the lattice. Based on symmetry [96] the effect of the field along [110] on the exchange constants is modelled through J 13 0 = J +ηz 0 = J 0 − δ(h) and J 24 0 = J −ηz 0 = J 0 + δ(h); for simplicity, we keep the other exchange constants J ij 0 unchanged (see Fig. 1b)). In order to detect the formation of a dense phase of monopoles in our Monte Carlo simulations we measure ρ and two more specific quantities: i− the average of the staggered O displacement along the z−axis, δu z stagg , that is sensitive to the O-displacement proposed by Sazonov and collaborators. We compute it as the average of δu z on up tetrahedra minus that on down tetrahedra (see Fig. 4a)); and ii− the order parameter, OP , for the double-layer crystal of single monopoles, calculated as the staggered charge on [100] planes made of up-tetrahedra. If we call Q up j the total charge in the j−th [100] plane of up tetrahedra, the OP is computed as: where L is the linear size of the system and we are counting two planes of up tetrahedra per unit cell. Fig. 5 shows the results obtained for the complete MeSI model of Eq. 7 as a function of temperature (filled symbols). We used a fixed field h/J ml = 13.4, with δ(h)/J ml = −0.5. In order to guarantee a spin ice phase at zero field we set J 0 /J ml = 1.1 > 1. The condition to destabilize the spin ice state in favour of a monopole phase at zero temperature is δ(h) < J ml − J 0 . With J 0 /J ml = 1.1, we make sure that a two-in-twoout state is favoured for h = 0 (δ(0) = 0), compatible with the observed Coulomb phase in Tb 2 Ti 2 O 7 . On the other hand, the value δ(h)/J ml = −0.5 ensures a single monopole phase at a finite field. We see the density of single monopoles saturating smoothly at ρ = 1 below T /J ml = 0.2. Since the intensity of the magnetic field was chosen in order that the α−spins would be saturated (and thus the magnetisation) for T /J ml < 1.4 this increase in ρ involves only the (antiferromagnetic) arrangements of β−spins. We can see that the staggered average of the O displacement along z−axis, δu z stagg , follows closely this behaviour, showing that the O-displacement along z−axis seem to coincide with that predicted in ref. 25. The negative value we needed to use for δ(h) is quite encouraging: it means that J ij = J +ηz in the link parallel to the field increases with field, while the one perpendicular (J −ηz , painted green in Fig. 4b)) decreases. The crystal contracts along the [110] field and expands in the direction perpendicular to it, in full accordance with the observed distortions under magnetic field [45,97].
In spite of the above, we notice that the specific heat C V (Fig. 5c), full symbols) shows only a broad Schottky anomaly on decreasing temperature, while the order parameter OP varies very little. This tells us that the spin ice-like ground state has changed into a dense monopole phase of incoherent β Q −chains, producing no spontaneous symmetry breaking. It is easy now to see that an effective interaction like the one in Eq. 6 with attraction between like charges (γ < 0) is the coupling needed to obtain the double monopole layer structure, since it favours the positive (negative) charges in contiguous β Q −chains to be next to each other (Fig. 4a)). It can be the result of second and third nearest neighbours exchange interactions [67], and may be related to the "half moons" in Tb 2 Ti 2 O 7 neutron scattering experiments [46,56]. Alternatively, the additional term can also be thought as an effective way to include the effect of the electric dipolar interactions, that have been proved to lead to the double-layered monopole crystal [54]. It is important to stress that their role here is not to stabilize single monopoles [54], but (more subtly) to favour a particular monopole arrangement. as a function of temperature and fixed field h|| [100], with h/J ml = 13.4 and L = 8. In order to guarantee a spin ice phase at zero field, J0/J ml = 1.1, δ(0)/J ml = 0, and δ(h)/J ml = −0.5. Curves with filled symbols where obtained using Eq. 7, while those with open symbols are the equivalent after adding an effective charge-charge interaction (Eq. 6) with γ/J ml = −0.08. a) Monopole density ρ (left) and average staggered O displacement along z−axis, δu z stagg (right). b) Order parameter. c) Specific heat. The inclusion of some attraction between like charges (Eq. 6) leads to an spontaneaus symetry breaking transition leading to the magnetic and structural phase shown in Fig. 4 .
The open symbol curves in Fig. 5 show the marked changes we measured after adding the monopole interac-tion term (Eq. 6) with γ/J ml = −0.08 to the extended MeSI model. We can see that ρ and δu z stagg reach saturation in a much sharper way. The abrupt jump in OP (reaching the value of 1), and the peak in C v (with an extra area under it) show that these sharp features are connected with the spontaneous symmetry breaking by the double monopole layer structure. In addition to the spin and monopole configurations, displayed on Fig. 4a), our model provides the lattice distortions linked to this magnetic structure (Fig. 4b)). As discussed before (Figs. 3b) and 4b)) and differently from Ref. 25, this displacement is not only vertical: the O ion tends to approach the triangular surface of each tetrahedron where the three spins point likewise (darkened in the figure), so as to reduce the value of the exchange constants along the corresponding links.
Given the big magnetic moments associated with Tb +3 , a brief consideration is needed regarding dipolar magnetic interactions. As discussed in Ref. 28, their effect will be twofold. Firstly, the preference of these interactions for two-in/two out states should be compensated by the huge magnetostriction of Tb 2 Ti 2 O 7 (i.e., the transition into a dense monopole phase would occur at higher fields/deformations than if no magnetic dipolar forces were included). Secondly, dipolar interactions would disfavour the proximity of like magnetic charges, demanding bigger values of |γ| (i.e., bigger next-nearest neighbours exchange interactions, or dipolar electric moments).

III. DISCUSSION
It is interesting to compare the resulting structures for some of the ground states which combine a maximum density of single monopoles and extensive residual entropy. We can now put together three pieces of information: the usual (spin) magnetic scattering, scattering from the (distorted) O −2 -lattice, and hypothetical scattering from magnetic charges. They were calculated from simulations at very low temperature, so that magnetic excitations are negligible and O −2 -ions are displaced only along the unit cell diagonals (see Section II C, and Supplementary Information II and III for details on the simulations and the precise definitions used in the Structure Factor calculations). Fig. 6 shows a comparison of the calculated structure factors within the [h, l, l] plane of reciprocal space. Three of the monopole phases (including the Polarizad Monopole Liquid -PML-studied in detail in reference 35) run along the rows of the table. The "scattering centres" (spins, magnetic charges, and the O −2 ions displaced from the centre of each tetrahedron) run along its columns. For the O −2 displacement we show only the diffuse part, removing the trivial contribution from the regular diamond lattice formed by the O −2 average position and k-dependent charge (see Supplementary Information III). On the other hand, Bragg scattering peaks due to static spins (polarised by the field or associated to the curl-free part of the magnetic moment in the crystal of single monopoles) are indicated schematically by full circles.
Monopole order progresses downwards in this table, as illustrated by the second column: broad maxima for the Monopole Liquid give place to pinch points in the PML and then to sharp Bragg peaks for the Zincblende Monopole Crystal. Regarding the Monopole Liquid, in spite of the maxima in the monopole channel, it shows no spin-spin correlations at all, which is also true for the O −2 displacements (first row in Fig. 6). The existence of these peaks in the charge channel may be counterintuitive given the total spin decorrelation. Monopolemonopole correlations in the ML are due to construction constraints, due to the underlying spins [24,29].
As previously mentioned in the text, the Zincblende Monopole Crystal shows pinch points both in the spin and the O −2 channel; strong Bragg peaks reflect the monopole correlations in the crystal. As the ML, the Polarized Monopole Liquid (middle row) has no spin or monopole long range order [35]. Notably, and differently from its unpolarized version, the PML has associated a gauge field that can be related either to spins, magnetic charges or displaced O −2 -ions. In principle, radiation interacting with any of these three particles could show the pinch points characterizing a Coulomb phase. The ability to detect effects related to the electric dipole on monopoles depends on its magnitude. While there are indications of the presence of such electric dipoles in Tb 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 [60,98], the O displacement δr has not been measured. Within our model, we can obtain it through J ml (Eq. 4), provided we know the coupling constantα and the elastic constant K. A rough estimation based on the experimental and numerical results obtained in Refs. [48,[99][100][101] gives a small J ml for Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 , on the order of 10 mK. In turn, this leads to δr ≈ 0.1 pm for these canonical spin ices. On the other hand, Jaubert and Moessner [54] estimate a bigger δr for Tb 2 Ti 2 O 7 (within the picometer range), similar to that observed in multiferroic materials. While this is still considerably small, new methods based on traditional XRD have been very recently proposed and used to observe distortions within this range in a strontium titanate oxyde [102]. The chances of a direct observation of the magnetoelastic phenomena we propose can increase if the efforts are first concentrated on compounds with a big coupling between magnetic and lattice degrees of freedom, starting with monopole crystals. The double monopole layer in Tb 2 Ti 2 O 7 could then be an excellent place to begin.
In summary, we have introduced an extension to the usual Hamiltonians used for studing Ising spin systems on pyrochlore oxides R 2 M 2 O 7 . The Magnetoelastic Spin Ice (MeSI) model includes the spin coupling to the lattice of central O −2 -ions in up and down tetrahedra, through the dependence of the superexchange constant J(δu) on the oxygen displacement (δu). We showed that, in the strong  [100] in the Polarized Monopole Liquid (PML), are indicated schematically by full circles. Three different gauge field types can be observed, associated to spins, magnetic charges, and oxygen displacements. While the Monopole liquid does not show pinch points in any case, the PML is most remarkable, reflecting the existence of a Coulomb phase in the three channels. The Fragmented Coulomb Spin Liquid [26,31,86] should reflect the existence of a Coulomb phase also through lattice distortions.
coupling limit, lattice distortions turn single monopoles (the excitations of the spin ice materials) into actual building blocks for novel ground states with maximum density of magnetic charges. Crucially, δu works as a dynamic, internal field; there is thus no explicit symmetry breaking, and all eight single monopoles are a priori equally probable in each tetrahedron.
This avenue to new ground states and novel physics is widened by an additional factor: the O −2 distortion implies an electric dipolar moment. This means that the distortions δu are not just "hidden" degrees of freedom that allow for the occurrence of new phases, but can be thought as probes to investigate the underlying magnetism, or, in a multiferroic fashion, to control the material.
We have presented some examples of the above. The first one is a Monopole Liquid ground state, stabilized for the first time with a Hamiltonian which physical bases. Including an attraction between magnetic monopoles of the same charge leads to "half moons" features in the spin structure factor of this liquid; this makes a direct link to the "spin slush" phases [30,61]. Remarkably, notwithstanding its simplicity, the MeSI model provides a unified framework that explains the zero field ground state measured in Tb 2 Ti 2 O 7 [55] and the double-layered monopole crystal at moderate fields [25]; this includes an improved version for the distortion on the O −2 -lattice proposed previously [25]. The classical treatment of distortions at low temperatures, albeit unrealistic, can be understood as a simple way to probe the magnetoelastic instabilities of this system. At zero field, the MeSI model recreates the phase diagram for single monopole spontaneous crystallization studied in [23,26] without resourcing to artificial constraints, and would allow to extend it to include double monopoles [24]. The spontaneous crystallization of a dense liquid of single monopoles into the Zincblende structure gave rise to the Fragmented Coulomb Spin Liquid [26,28,31]. As stressed before, there is a close parallelism between some electric and magnetic phenomena in frustrated Ising pyrochlores [53]. Our access to the elastic degrees of freedom provides another layer of complexity to this, by showing that the O −2 displacement δu in the FCSL phase is related (like the spin) to a gauge field. Perhaps more singular is the case of the Polarised Monopole Liquid [35] (i.e., the monopole liquid with an applied magnetic field along [100]). This disordered state is a Coulomb phase from the point of view of three different degrees of freedom: spins, magnetic monopoles and elastic distortions. As with the FCSL, pinch points could be detected using diffuse neutron scattering or simply xray or electron diffraction.
Although here we have concentrated mainly on the magnetic degrees of freedom and on the strong coupling limit, the MeSI model opens perspectives of research in other grounds. For instance, the weak coupling regime can be used to describe in a combined way (spins and lattice distortions) some of the physics of spin ice materials [53,57,58,60]: their true ground state [48,49,[103][104][105][106][107], the effect of uniaxial pressure [48,108], and the new phases that appear under applied field [49,74,83,84,109].

IV. ACKNOWLEDGEMENTS
This work was supported by Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT) through grants PICT 2013-2004, PICT 2014-2618 and PICT 2017-2347, and Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) through grant PIP 0446. Part of this project was carried out within the framework of a Max-Planck independent research group on strongly correlated systems.

A. Derivation of the Hamiltonian
Most of the results shown in this work have been derived from Eq. 2 of the main text, obtained after rewriting the magnetoelastic exchange energy for small O −2 distortions δu under certain assumptions. Here we will make explicit this procedure. The energy we are concerned with corresponds to the second term in Eq. 1 of the main text; for a single tetrahedron it can be written as 1 Note that the {S i } are pseudospins; a positive value of the superexchange constant J ij , for instance, would favour pairs of pseudospins of different sign, which translates into a ferromagnetic-like arrangement in the tetrahedron (one spin in, the other pointing out) [9]. As justified in the text, we assume that the superexchange interaction between the nearest-neighbour spins S i and S j is mediated by the O −2 -ion near the centre of the spin arrangement, and that it is a function of the angle θ ij formed by the position vectors of the spins relative to this ion, J ij = J(θ ij ) (Fig. 7). We will now obtain a functional dependence of the couplings J ij on the displacement of the oxygen δr, for very small displacements (|δu| ≡ |δr|/r nn 1); we assume, for simplicity, that the spins remain in fixed positions of the lattice.
We concentrate first on an up tetrahedron. When the oxygen is at the central position, the θ ij angles are all equal to θ = arccos(− 1 3 ) = 109.47°. To lowest order, the cosine is linear in δθ when expanded around this value. Let us, for instance consider the angle θ 13 between the relative position vectors of S 1 and S 3 with respect to the displaced oxygen (see Fig. 7). The positions involved (measured from the exact centre of the tetrahedron, in units of r nn ) are u 1 = √ 2/4 × (1, 1, 1) for spin 1; u 3 = √ 2/4×(−1, −1, 1) for spin 3, and δu = (δu x , δu y , δu z ) for the O −2 -ion. The expression for θ 13 can be obtained via the scalar product between its defining vectors: . From this we get: .
Taylor expanding this expression to first order in δu x , δu y , and δu z , it is easy to see that where C = 8 3/2 /9. In this first order approximation only the displacement along z -which is the only coordinate that has the same value for spins 1 and 3-remains. It may be useful to remember that if this displacement has a positive value, it yields a negative contribution to Eq. 9, implying a more obtuse angle. Following Goodenough-Kanamori-Anderson rules for superexchange, we can then expect that the O −2 -ion displacement decreases the value of the exchange constant of the magnetic link it approaches (it makes it less ferromagnetic, or more antiferromagnetic). This is illustrated in Figs. 1b), 3 and 4b) of the main text, where we can see that pairs of spins closer to the O −2 -ion (linked by green segments) tend to point all out or all in.
Proceeding analogously for the angle between spins S 2 and S 4 , the common coordinate is again z, but this time with a negative sign; it yields a term with a positive sign: This can be repeated, for the other angles obtaining: For each case, only the coordinate with the same value for both spins considered produces a first order term, and this term has opposite sign with respect to the value of the coordinate. The same equations but with opposite sign in the second term would be obtained for a down tetrahedron (for instance, a positive δu z brings the O −2 nearer the link connecting spins S 1 and S 3 in an up tetrahedron (Fig. 7), but it does the opposite in a down one).
The previous calculation was purely based on geometry. Resuming physical grounds, we now expand assumed positive [53,54]. J 0 is the zeroth order approximation, J 0 = J ij (δu = 0), shared by all the links in the absence of structural symmetry breaking (see for instance the case of a field applied along [110], Eq.5 in Section V of the main text). If, for instance, we sum the Hamiltonian terms for an up tetrahedron involving δu x , we obtain −αδu x (S 1 S 2 − S 3 S 4 ), which is just the first addend of the scalar product in the second term of Eq. 2 of the main text. The addends involving J 0 have been gathered in the third term in this Hamiltonian, which amounts to the usual nearest neighbours spin Hamiltonian [13]. The first term in Eq. 2 of the main text accounts for the elastic energy; in the absence of magnetoelastic coupling, this term ensures that the equilibrium position of the central oxygen is given by δu = (0, 0, 0) (i.e., at the centre of the tetrahedron).
B. Details on the numerical simulations.
We provide here some details on the simulations used in the main text to study the equilibrium properties of the Ising pyrochlores. We simulated L × L × L conventional cubic cells of the pyrochlore lattice (16 × L × L × L spins) with Metropolis and Conserved Monopoles [23] algorithms. In all cases the boundary conditions were set to be periodic along the cubic primitive vectors. The specific heat and susceptibilities were calculated using the fluctuations of the corresponding quantities.

Metropolis algorithm
Given that the MeSI model is a composite system, we separated each step of the simulation into a single spin flip or single elastic O −2 movement. The Metropolis algorithm was then implemented, using Eq. 2 of the main text to evaluate the probability to accept or reject a given spin flip or O −2 -ion movement. Typical sizes were L = 8 (N = 8192 spins, and N t = 2048 O −2 -ions). To simulate the elastic distortions we considered displacements of the O −2 -ions using spherical coordinates, with δu chosen randomly in a distribution from 0 to a temperature dependent maximum δu max (T ). The latter is introduced so as to make more efficient the algorithm [42]. We considered the relaxation times of the magnetic degrees of freedom to be much shorter than the corresponding elastic ones, thus, within each elastic move we performed a complete magnetic Monte Carlo step. We verified that this election has no effect on the equilibrium properties of the system. After reaching equilibrium, we averaged the data over 50 independent runs, taking 2 × 10 4 time-steps at each field and temperature.

Conserved Monopoles Algorithm
Energy minimisation in a single monopole site requires the O −2 -ion to be situated along the 111 diagonals, approaching the tetrahedron face with all spin pairs in antiferromagnetic-like configurations (see Fig. 3b) of the main text). This fact makes easy, in a ground state, to infer the O −2 distortions {δu} if the magnetic degrees of freedom {S i } are known. We have profited from this fact to calculate the structure factors for different dense single monopole ground states with better statistics and in larger lattices with a minimum computational effort. Using the conserved monopoles algorithm [23,74] we travelled along the desired magnetic ensamble, and then deduced the elastic coordinates for each spin configuration. While each configuration had in truth two neutral sites, the density of monopoles is very near ρ = 1 for the system sizes used (for L = 8 the proportion between neutral sites and single monopoles is 2/4094).

C. Calculation of the different structure factors
The simulated neutron structure factors have been calculated following the expression: where i and j sweep the pyrochlore lattice, N is the number of spins, and ... represents thermal average (in this case, that of the product of pseudospins at sites i, j).
The spin quantization directions are given by {μ i } (parallel to the 111 directions). Then, µ ⊥ i is the component ofμ i of the spin S i = S iμi at site i perpendicular to the scattering wave vector k: As in Refs. 29 and 35, we also calculated a structure factor associated with magnetic charges, through the Fourier transform for the charge-charge correlation function: where the greek indices α, β now sweep the sites of the diamond lattice where the monopoles live, N/2 is the number of tetrahedra, Q α represents the topological charge at position r α , and r αβ is the distance between monopoles. Finally, we have measured the diffuse structure factor associated to the displaced O −2 -ions near the ground state of the system. Assuming an atomic form factor unity, we calculated: where q av (k) = e ik·δrα is an average, k−dependent O −2 "charge" in the perfect diamond lattice. These functions have been obtained after thermal averages over sets composed of 500−1000 independent configurations for a system size L = 8.