Monopole matter from magnetoelastic coupling in the Ising pyrochlore

Ising models on a pyrochlore oxide lattice have become 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 the double-layered monopole crystal observed in Tb2Ti2O7. Additionally, the dipolar electric moment of single monopoles emerges as a probe for magnetism. As an example we show that some Coulomb phases could, in principle, be detected through pinch points associated with O−2-ion displacements. Pyrochlore oxide lattices are home to many geometrically frustrated magnetic systems, including the celebrated spin ice. The authors present an Ising magnetic model on this structure, with a magnetoelastic term coupling the spin lattice to the O−2 ions which mediate superexchange. This model explains the presence of exotic forms of order found in previous experiments, allows stabilisation of some long-sought phases, and signals lattice distortions as a new probe for the complex magnetism of these materials.

I t is remarkable that the Ising model, one of the simplest interacting systems in condensed matter physics, can lead to phenomena in geometrically frustrated magnetism that have kept researchers interested for the past decades. The strategic choice of the spin lattice structure (such that pairwise interactions compete rather than collaborate) is 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][2][3][4][5][6][7] , among other effects.
The pyrochlore structure (Fig. 1a) is a prominent instance of 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 non-negligible 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 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 spin configurations related to single monopoles in spin ice 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 stabilisation in pyrochlore oxides. We begin by introducing the Magnetoelastic Spin Ice (MeSI) model, by modifying the nearest neighbours spin ice Hamiltonian in a simple way to include a coupling to the lattice of O −2 -ions sitting near the centre of tetrahedra (see Fig. 1). In the strong coupling regime, lattice distortions turn the eight types of single monopoles into stable, atomic-like constituents of novel ground states. We then show how the MeSI Hamiltonian stabilises a Monopole Liquid; this massively degenerate perfect paramagnet is 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,55 , 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 . 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 this way, we will see that in the Zincblende Monopole Crystal the deformed O lattice fluctuates with the fragmented magnetic moments. Also, in spite of its simplicity, and building on the previous works of Jaubert and Moessner 54 and Sazonov and collaborators 25 , our 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 56 , and gives some clues on the half-moons measured in neutron diffraction patterns at finite energy 46,57 .
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 electrical In the oxide pyrochlore structure the Ising spins (coloured arrows) occupy the shared vertices of tetrahedra, pointing either towards or away from their centre. Like in an electrostatic field, they can be associated with a flux density. Each spin configuration can be related to a magnetic charge at the centre of the tetrahedron, through the net flux into it. For instance, the three-in/one-out configuration in the "up" tetrahedron embedded in a cube has a positive single monopole (small green sphere), while the one-in/three-out in a "down'' tetrahedron (coloured blue) has a negative single monopole (small red sphere). Neutral tetrahedra (with no spheres) are characterised by twoin/two-out configurations. The exponentially degenerate ground state of Spin Ice is built from these magnetically neutral blocks. Double charges, with all spins pointing in or all-out, are also possible: a negative double charge is represented here by a big red sphere. Unfrustrated crystals of double charges alternating in sign are known ground states of various materials. Single monopoles are the low energy excitations of both spin ice and of these crystals. However, we have found that a coupling between elastic and magnetic degrees of freedom can stabilise a ground state of single monopoles that preserves all the symmetries of the lattice they sit in (i.e., a "liquid" of monopoles). b The oxygen ion (small cyan sphere) that mediates superexchange between nearest neighbours spins (generally associated to magnetic rare earth ions, R, purple spheres) minimises the elastic energy at the centre of the tetrahedron. Its displacement δr toward the + z − link decreases the exchange constant value between spins 1 and 3 (connected by a green line); the exchange constant in the opposite link (−z − link, red line) is increased by this distortion, and the other four (black lines) J values remain unchanged to first order. δr thus favours certain spin configurations, while a given configuration conditions the distortions of the oxygen lattice.
properties of Ising pyrochlores, or teach us how to probe other properties through them (as it has been done in some pioneering works in spin-ice [58][59][60][61] ).

Results
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-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); these monopoles are represented 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-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 neighbour magnetic interactions, our model Hamiltonian can be written as Here δu ≡ δr/r nn , with r nn the nearest neighbour 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 neighbours superexchange energy associated to each pair. It can also be labelled 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 Note 1.
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 Note 1, to first order J ±m (δu) is only affected by the mcomponent of δu: J ± ηm ðδuÞ ¼ J 0 Àð ± Þηαδu m , where η takes the value + 1 (− 1) for up (down) tetrahedra. The constantα ∂J ± m ∂δu m δu¼0 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 Note 1) 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 nearest-neighbour 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 "Doublelayered crystal of single monopoles") the latter would be replaced by a sum involving the exchange constant matrix J ij 0 : The middle term in the MeSI model is central to this work, as it contains the (linearised) 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.
Stabilisation of a dense fluid of single monopoles: the monopole liquid. Models of interacting entities, even simple ones, are seldomly exactly solvable. It is then a surprise that the threedimensional 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) stabilises 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 stabilise 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. Figure 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 (without counting 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 an increase in 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 additional 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 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-neighbour spin ice term proportional to J 0 and the four-spin 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 nearest 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, ð6Þ where = 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  2), with coupled elastic and magnetic degrees of freedom (symbols). This model supports a liquid of single magnetic monopoles as a ground state. In order to see this, the exchange constant in the absence of distortions was set to J 0 = 0, the magnetic field to h = 0, and the energy scale associated with the stability of magnetic monopoles to J ml ¼ 3α 2 =K ¼ 1K. The cubic lattice we simulated had N= 8192 spins and N/2 moving O −2 ions. By integrating out the lattice degrees of freedom, this compound system can be transformed into an effective magnetic model (Eq. (5)) studied previously 29 . The analytical results from Barry and Wu 66 for the effective model are also plotted (full lines). a Left axis: density of monopoles ρ as a function of temperature. ρ saturates at 1 single monopole per tetrahedron for T << J ml , while the system retains all the symmetries it had at high temperature. Right axis: inverse of the magnetic susceptibility χ vs. temperature. 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 . b Specific heat per unit spin C V (left) and the mean quadratic deviation δu 2 (right axis) vs temperature. Both curves make evident that we are dealing with a coupled composite system, with magnetic excitations and lattice distortions. 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. Also, δu 2 (T) is not just a straight line: the O −2 equilibrium position depends on the spin configuration in δO. 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 55 , with single monopole density ρ ≈ 0.35 and attraction between like-monopoles. We have calculated the neutron structure factor I Spin (k) for the ML in the presence of nearest neighbour attraction between like charges as per Eq. (6) with γ < 0, ρ = 1, T << |γ| < J ml (see Supplementary Note 3 for details). The diffuse pattern we obtain can be understood as the result of merging the different half-moons observed in refs. 30 and 55 , with their features widening due to a higher density. While half-moons are usually detected at finite energy 57,68 , the ML with like-attraction is a new instance (together with refs. 30,55,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 "Doublelayered crystal of single monopoles"), a compound which also shows half-moons in its neutron structure factor 46,57 .
The fragmented Coulomb spin liquid (FCSL): correlated magnetic and dipolar electric fluctuations. 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 fields [33][34][35]76,83 , internal fields 26,28,31,84 , 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 , characteristic of a Coulomb phase 4 . A number of FCSL have been recently achieved experimentally in a pyrochlore lattice 31,84 . 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 "Stabilisation of a dense fluid of single monopoles: the Monopole Liquid" into a ZnMC when the temperature is lowered. 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,85 . Figure 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 magnetic charges in a lattice 23 , analogous to the one obtained for electric charges in a lattice 86 , 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 stabilised 24,87 . Defining the Fig. 3 Zincblende monopole crystal (J 0 = 0 and γ/J ml = 0.2). An attraction between monopoles of different signs can make the Monopole Liquid crystallise in the Zincblende structure, with spontaneous symmetry breaking. The ordered charges in the crystal contrast with the partial disorder of the spins that produce this crystal structure 23 , in what is known as the magnetic moment fragmentation 26 . a Specific heat (C V , right axis) and density of single monopoles (ρ, left axis) across the temperature at which the crystallisation transition occurs. b The order parameter for crystal formation (the staggered monopole density ρ S , left axis) and its fluctuations χ S (right) measured in the same temperature range. We have subtracted the contribution of the pure vibrational degrees of freedom from C V . 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 links (with the three spins on the corners sharing the same sign), coloured green here. The fluctuation of the fragmented magnetic moment implies that δr is also fluctuating. 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 stabilised 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 hence the electric dipole moment, remain 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 "Discussion", where we also summarise 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.
Double-layered crystal of single monopoles. Among the complex physics of Tb 2 Ti 2 O 7 there is a clear experimental fact: upon applying an external field h parallel to the [110] crystallographic direction, an order of alternate double layers of positive and negative monopoles is induced perpendicular to [001] 25,88 (see Fig. 4a). To justify this charge order, Jaubert and Moessner 54 explored a classical model. The mechanism involved the longrange 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 "all-in/all-out" phase into the bi-layer when applying a [110] magnetic field.
The MeSI model constitutes an alternative to this intrinsic mechanism, the first ever proposed that could stabilise a monopole phase 54 . While it obviously cannot take into account all the complexity observed in Tb 2 Ti 2 O 7 56,64,89-92 , it is an improvement over the previous proposal, since now both magnetic and elastic degrees of freedom are considered on an equal footing. Furthermore, the model provides a unified explanation for the ground state observed at zero field (a Coulomb phase 56,93 ), 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,57 .
The application of a strong magnetic field along [110] does not fully order the Ising pyrochlore lattice. Spins on α-chainsrunning along [110], represented by purple arrows in Fig. 4aare completely polarised 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 77 ; 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 issue of charge coherence, the non-trivial question we should answer concerns magnetic charge stability. Why would a sufficiently strong field h||[110] change the ground state of Tb 2 Ti 2 O 7 from a subset of the 2-in 2-out manifold to that of a dense phase of single monopoles 25 ? 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 The spin lattice can be divided in: i-α − spins (painted purple) polarised by h; ii-β − spins (yellow) perpendicular to h and building up antiferromagnetic β − chains (in-in/out-out). Monopoles linked by β-spins alternate in sign, forming β Q − chains of charge. α − spins decouple consecutive β Q − chains but, as proposed here, an effective nearest-neighbours attraction between like monopoles (Eq. (6) with γ < 0) can favour order such that the charge in two tetrahedra linked by an α − spin is the same. The main question our model is able to answer concerns the nature of the force that stabilises the single monopoles occupying each tetrahedron when the field increases from zero. b Structural configuration for the O −2 -ion near the centre of each tetrahedra. Through the magnetoelastic coupling, monopoles are stabilised by the displacement of the central oxygen (small cyan spheres) along 〈111〉 directions. These distortions (which should be compared with those proposed on Fig. 6 in ref. 25 decrease the exchange energy along the three bonds of the triangular face (painted green) approached by the O −2 ions, favouring "three-in", or "three-out" configurations in these triangles (see a)). The spins form "one-in/one-out" configurations on the other bonds (painted red), where the exchange constant increases. The exchange constants along ±z (painted in slightly darker colours) are further modified by magnetostriction, which triggers the O −2 displacement. 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 that is big enough to overcome the energy associated to the usual nearest neighbours term, proportional to J 0 . Alternatively we will now investigate the effect of h on the lattice structure, and on the spin lattice mediated by the magnetoelastic couplingα.
We are not the first to notice the possible importance of the giant magnetostriction observed for h||[110] 45 for stabilising 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 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 94 the effect of the field along [110] on the exchange constants is modelled through 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, hδu z stagg i, that is sensitive to the O-displacement proposed by Sazonov and collaborators, computed 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. Figure 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 destabilise 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-two-out 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. The density of single monopoles saturates 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, hδu z stagg i, 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,95 .  7)) for a fixed magnetic field h||[110], with h/J ml = 13.4 and 8192 spins. The extended model indirectly takes into account magnetostriction effects using inhomogeneous exchange constants J ij 0 , modifying J 0 according to symmetry: Fig. 1b)); the other exchange constants for the undisplaced O −2 were kept constant at the value J 0 . In order to guarantee a spin ice-like phase at zero field, J 0 /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 chargecharge interaction (Eq. (6)) with γ/J ml =− 0.08. a Monopole density ρ (left) and average staggered O displacement along z − axis, hδu z stagg i (right axis) as a function of temperature. b Order parameter for the crystal (OP) as a function of temperature. c Specific heat (C V ) vs temperature. The attraction between like charges (Eq. (6)) leads to an spontaneous symmetry breaking transition into the magnetic and structural phase shown in Fig. 4. It is reflected in the abrupt changes seen in all the curves near T/J ml ≈ 0.5.
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 charges of equal sign 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,57 . 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 doublelayered monopole crystal 54 . It is important to stress that their role here is not to stabilise single monopoles 54 , but (more subtly) to favour a particular monopole arrangement.
The open symbol curves in Fig. 5 show the marked changes we measured after adding the monopole interaction term (Eq. (6)) with γ/J ml = −0.08 to the extended MeSI model. We can see that ρ and hδu z stagg i 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).

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 "The Fragmented Coulomb Spin Liquid (FCSL): correlated magnetic and dipolar electric uctuations", and Supplementary Notes 2 and 3 for details on the simulations and the precise definitions used in the Structure Factor calculations). Figure 6 shows a comparison of the calculated structure factors within the [h, l, l] plane of reciprocal space, laid out forming an array. Three of the monopole phases (including the Polarised Monopole Liquid -PMLstudied in detail in a previous work 35 ) run along the rows. 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 Note 3). 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 figure array, 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 Polarised Monopole Liquid (middle row) has no spin or monopole long-range order 35 . Notably, and differently from its unpolarised version, the PML has a gauge field associated that can be related to either spins, magnetic charges or displaced O −2 -ions. In principle, radiation interacting with any of these three particles could show the pinch points characterising 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 61,96 , the O displacement δr has not been measured experimentally. 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,[97][98][99] 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 picometre 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 oxide 100 . 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 starting point.
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 show that, in the strong coupling limit, lattice distortions turn single monopoles (the excitations of the spin ice materials) into actual building blocks of 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, stabilised for the first time with a Hamiltonian with 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,55 . 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 56 and the double-layered monopole crystal at moderate fields 25 . This includes an improved version for the previously proposed distortion of the O −2 -lattice 25 . The classical treatment of distortions at low temperatures presented here, 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 crystallisation studied in 23,26 without resourcing to artificial constraints, and can be suitably extended to include double monopoles 24 . The spontaneous crystallisation 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 . Access to the elastic degrees of freedom provides another layer of complexity, by showing that the O −2 displacement δu in the FCSL phase is, like spins, related 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 by means of x-ray or electron diffraction.
Although we have concentrated our discussion mainly on the magnetic degrees of freedom and on the strong coupling limit, the MeSI model opens perspectives of research in other grounds. For The structure factor shows this in a very illustrative way. In this array, the nature of the scattering centre (spins, single monopoles, or the oxygen ions near the centres of the tetrahedra) varies along the horizontal axis (see Supplementary Note 3 for details on the structure factor calculation). The names of three of the studied phases are indicated on the left, with increasing charge order (while not necessarily spin order) progressing downwards. For scattering from O −2 ions we have omitted the Bragg peaks corresponding to their position in the regular pyrochlore lattice; on the other hand, Bragg peaks from the fragmented divergence-full component of the magnetic moment in the Zincblende Monopole Crystal (ZnMC), or by the spins aligned by the applied magnetic field along [100] in the Polarised Monopole Liquid (PML), are indicated schematically by full circles. Notably, three different types of gauge field can be observed, associated with 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. In principle, the Fragmented Coulomb Spin Liquid 26,31,84 should reflect the existence of a Coulomb phase not only through the magnetic moments, but also through lattice distortions.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The codes used to generate the results presented in this manuscript are not widely available. However, any parties wishing to implement, wholly or partly, the framework presented in this manuscript are invited to contact the authors, who will make every effort to assist with implementation of this model where reasonably possible.