# Spin ice under pressure: symmetry enhancement and infinite order multicriticality

Ludovic D.C. Jaubert, J.T. Chalker, Peter C.W. Holdsworth, R. Moessner

► To cite this version:

Ludovic D.C. Jaubert, J.T. Chalker, Peter C.W. Holdsworth, R. Moessner. Spin ice under pressure: symmetry enhancement and infinite order multicriticality. Physical Review Letters, 2010, 8 (105), pp.087201. 10.1103/PhysRevLett.105.087201. ensl-00462955

HAL Id: ensl-00462955

https://ens-lyon.hal.science/ensl-00462955

Submitted on 10 Mar 2010

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Spin ice under pressure: symmetry enhancement and infinite order multicriticality

L. D. C. Jaubert,¹ J. T. Chalker,² P. C. W. Holdsworth,³ and R. Moessner¹

¹Max-Planck-Institut für Physik komplexer Systeme, 01187 Dresden, Germany.

²Theoretical Physics, Oxford University, 1 Keble Road, Oxford, OX1 3NP, United Kingdom. and

³Laboratoire de Physique, École Normale Supérieure de Lyon,
Université de Lyon, CNRS, 46 Allée d'Italie, 69364 Lyon Cedex 07, France.

(Dated: March 10, 2010)

We study the low-temperature behaviour of spin ice when uniaxial pressure induces a tetragonal distortion. There is a phase transition between a Coulomb liquid and a fully magnetised phase. Unusually, it combines features of discontinuous and continuous transitions: the order parameter exhibits a jump, but this is accompanied by a divergent susceptibility and vanishing domain wall tension. All these aspects can be understood as a consequence of an emergent SU(2) symmetry at the critical point. We map out a possible experimental realisation.

PACS numbers: 05.50.+q 05.70.Jk 75.10.Hk 75.40.Cx 75.40.Mg

One fascinating aspect of condensed matter physics is the extent to which the nature and symmetries of emergent low-energy degrees of freedom can be independent of the high energy ones from which they derive. While many types of low energy degrees of freedom are possible, systems in which these correspond to a gauge field are still rare.

In this paper we study this theme in a setting which brings forth an unusual sequence of low energy fluctuations and symmetries: we consider spin ice, a frustrated Ising magnet on the pyrochlore lattice, whose low temperature behaviour is well-described by a gauge field, an emergent magnetostatics, representing a so-called Coulomb phase where correlations decay algebraically as $r^{-3}$. Transitions out of such a phase have attracted a great deal of attention recently [1–4].

We ask what happens when the host crystal is subjected to uniaxial pressure. This is a powerful probe of correlated matter, which has previously produced interesting information in frustrated magnets closely related to spin ice [5, 6]. Specifically, we consider the case of pressure-induced strain that lowers the crystal symmetry from cubic to tetragonal. We find that this induces a highly unusual symmetry-breaking phase transition out of the Coulomb phase: it is characterised by discontinuities in quantities such as the magnetisation, without being first order. Rather, it occurs at an infinite order multicritical point and is accompanied by a divergent susceptibility. We demonstrate both these features using Monte Carlo simulations. Further salient properties of the critical point include the absence of a domain wall tension, along with spin correlations that vanish in a plane perpendicular to the strain-induced tetragonal axis. These should be observable experimentally, and we discuss the scope for realising such a transition in the laboratory.

The unusual nature of the transition is a consequence of an emergent SU(2) symmetry: as the pressure reduces the symmetry of the external space, an enhanced internal symmetry appears. We show this by means of an exact

solution in three dimensions at the critical temperature and a mapping to a quantum phase transition in 2 + 1 dimensions, between Ising and XY anisotropy in a spin-1/2 ferromagnet.

Transitions with analogous features have been studied previously in the context of ferroelectrics [7], using versions of the model for potassium dihydrogen phosphate (KDP) introduced by Slater [8], which in its two dimensional form is equivalent to the 6-vertex model and exactly solved [9]. Spin ice has the potential to provide a much cleaner realisation of this physics than the transition in ferroelectrics, since its magnetic degrees of freedom are much more accurately represented by Ising variables supported on a rigid lattice.

Spin ice is well modelled by Ising spins $S_i$ parallel to their local easy-axis on the pyrochlore lattice [10]. At sufficiently low temperature we can limit ourselves to the highly degenerate ground state ensemble with two spins pointing in and two out of each tetrahedron (ice-rules) and consider effective nearest neighbour interactions [11, 12]. After coarse-graining the magnetisation, the discrete ice-rules can be written in terms of a continuous divergence free condition from which emerge the algebraic correlations [13, 14]. Our goal is to lift the degeneracy of the resulting Coulomb phase via the exchange modulation

$$ \mathcal{H} = - \sum_{\langle i,j \rangle} J_{ij} S_i \cdot S_j \quad (1) $$

where $J_{ij} = J - \delta > 0$ for bonds on (001) planes perpendicular to the strain axis, and $J_{ij} = J$ otherwise (see Fig. 1). For $\delta > 0$, the $\mathbb{Z}_2$ symmetry of this Hamiltonian is spontaneously broken at temperature $T_c \sim \delta$ in favour of a state magnetised parallel to the [001] direction. If the strain were very large, so that $T_c \sim J$, ordering would be from a conventional paramagnetic state and in the standard three-dimensional Ising universality class. In contrast, for realistic, small strains $T_c$ is much less than $J$ and ordering is from the Coulomb phase, with
the striking features we describe.

<figure><img src="image_2.png" /><figcaption>FIG. 1: A tetrahedron of the pyrochlore lattice with exchange $J$ and $J-\delta$ as indicated: the dashed bonds lie in (001) planes. One of the two ground states for $\delta > 0$ is shown.</figcaption></figure>

We first present results from Monte Carlo simulations, done with a worm algorithm[15, 16] in which loop updates preserve the divergence free nature of the magnetisation. The temperature dependence of the magnetisation $M$ and of the inverse susceptibility $\chi^{-1}$ are shown in Fig. 2, upper panel: remarkably, $M$ takes its saturation value $M_{sat}$ at all temperatures below $T_c$, but $\chi$ diverges as $T_c$ is approached from above. To better characterise the transition, we also examine the probability distribution function (PDF) of the order parameter, shown in Fig. 2, lower panel. This is a standard diagnostic: two peaks are expected in the PDF at the transition point when this is first order, but only one peak if it is continuous [17]. Again, we find unconventional behaviour: a PDF that is uniform over all values of $M$ at $T_c$.

An exactly uniform order parameter distribution arises within Landau theory as a limiting case. The free energy density $G(M)$ close to a multicritical point of order $n$ is

$$G(M) = \frac{\alpha_2}{2}(T - T_c)M^2 + \sum_{n'=n}^{+\infty} \frac{\alpha_{2n'}}{2n'}M^{2n'} \quad (2)$$

In the limit $n \to \infty$, it is independent of $M$ at $T = T_c$. Moreover, the order parameter exponent takes the value $\beta = 1/(2n-2)$, generating a jump for $n \to \infty$. Also, the susceptibility ($\chi \sim A_{\pm}|(T - T_c)/T_c|^{-\gamma}$ for $T \ge T_c$) has an exponent $\gamma = 1$ independent of $n$ and an amplitude ratio $A_{-}/A_{+} = 1/(2n-2)$ that vanishes as $n \to \infty$. All these are borne out by our simulations. This success of mean field theory reflects the exceptional value of the upper critical dimension, $d_{c2} = 1$, which follows from standard arguments with allowance for the anisotropic spatial scaling, see Eq 5 below.

In previous discussions of the ferroelectric KDP model [7, 8, 18–21] some distinctive aspects have been recognised. Indeed, Slater’s original approximation [8] yields a free energy independent of the order parameter at the transition, which on this basis has been called infinite order [22]. Such a conclusion is evidently very delicate [18], and corrections to the approximation or physical perturbations have the potential to convert the transition into a conventional one: either first order or continuous. This is in fact the case for the material KDP

itself, which has a first order transition that can be driven through a tricritical point under pressure [7]. While in two dimensions, an exact solution of the KDP model corroborates the main features of Slater’s results [9, 21, 23], it is not a priori clear – not least in view of the famous difficulties in generalising Onsager’s solution of the Ising model – what this tells us for higher dimension. Here we present simulations and exact results for the transition in three dimensions. In addition, we show how to realise it in a magnetic material.

<figure><img src="image_3.png" /><figcaption>FIG. 2: (Color online) Top: Simulation data for $M$ (blue circles) and $(\chi T_c)^{-1}$ (red triangles) vs. $T/T_c$. Bottom: PDF of $M/M_{sat}$ at different temperatures. For $T < T_c$ (red) the PDF has two sharp peaks at $M = \pm M_{sat}$. For $T = T_c$ (black) the PDF is broad in small samples (dashed black line for $L_{\perp} = L_z = 4$) and almost flat in larger samples (solid black line for $L_{\perp} = 4$ and $L_z = 18$). For $T > T_c$ (blue) the PDF is Gaussian. (Values at $T \ne T_c$ scaled to fit on vertical axis.)</figcaption></figure>

The ice rules impose the same value of the magnetisation $M$ in all (001) planes: any configuration can be mapped onto an ensemble of strings of down spins spanning the system from top to bottom, with the convention that the string vacuum is the state in which all spins are up [18]. The number of strings defines (i) the total magnetisation and (ii) a topological sector, so that the former can be used to label the latter. The equiprobability of all sectors (see Fig. 2) already suggests the absence of interactions between strings. This observation can be made precise using the transfer matrix that acts between two adjacent (001) layers of the lattice [24]. Since this transfer matrix is a direct product of factors $T$ representing separate tetrahedra, we first discuss a single tetrahedron. A single string can enter the tetrahedron at either site in one (001) layer and leave it at either site in the other layer, with energy cost $4\delta/3$. This sector is therefore rep-
resented in **T** by a 2 × 2 block in which all entries are the Boltzmann factor Γ = e⁻⁴δβ/³. By contrast, zero or two strings in a tetrahedron cost no energy and impose the next configuration locally, generating two 1 × 1 blocks with unit entries in **T**, which hence has the form

$$ \mathbf{T} = \begin{pmatrix} 1 & & & \\ & \Gamma & \Gamma & \\ & \Gamma & \Gamma & \\ & & 1 & \end{pmatrix}. \qquad (3) $$

The maximum eigenvalues of **T** lie for $T < T_c = \frac{4\delta}{3\ln 2}$ in the two fully magnetised sectors, but for $T = T_c$ all sectors have the same maximal eigenvalue. In addition, the associated eigenvector in the one-string sector is (1, 1). Explicit calculation of the transfer matrix for a complete lattice is complicated by the fact that the repeat unit in the [001] direction involves four different layers. Nevertheless, its maximal eigenvalues inherit precisely the properties described for **T**, and at $T = T_c$ the associated eigenvectors in each sector give equal weight to all arrangements of strings on a (001) layer. The maximal eigenvalues determine the physical properties for a system of size $L_z \times L_\perp^2$ if the [001] dimension satisfies $L_z \gg L_\perp^2$. As a result: (i) below $T_c$ the magnetisation is saturated; and (ii) at $T_c$ all sectors are equiprobable and all configurations within a sector are equiprobable, reflecting the absence of interactions between strings [16]. By contrast, for a system with $L_z \sim L_\perp$ subleading eigenvalues also contribute and the PDF for *M* is rounded, as is apparent in Fig. 2.

The transfer matrix for the complete lattice can alternatively be thought of as the evolution operator in imaginary time for a two-dimensional quantum XXZ model (in analogy with Eq. (8) of [9]) with Hamiltonian

$$ \mathcal{H}_Q = -J \sum_{\langle i,j \rangle} (s_i^x s_j^x + s_i^y s_j^y + \Delta s_i^z s_j^z). \qquad (4) $$

The strings thus denote the world lines of $s_i^z = 1/2$ spins. In this language the decomposition into sectors is a consequence of the conservation of total $s^z$, and is due to a U(1) symmetry. The ground state of $\mathcal{H}$ has ferromagnetic order of the quantum spins for all $\Delta$, with an orientation for the magnetisation that depends on $\Delta$. For $\Delta < 1$, spins lie in the x-y plane and the state has zero-point fluctuations, representing the Coulomb phase. For $\Delta > 1$, spins are aligned along the z-axis and the state has no fluctuations, representing the low temperature phase of the classical system. The equiprobability of string sectors at the critical point corresponds here to the degeneracy of the $N + 1$ ground states of the isotropic Heisenberg ferromagnet with $N$ spins, for all values of the total magnetisation $\Sigma^z = \{-N/2, \dots, N/2\}$. We have confirmed that this isotropy is a property of the full transfer matrix, and not only of the leading eigenvectors, by checking that the matrix commutes with the total spin raising and lowering operators, $\Sigma^\pm = \sum_j (s_j^x \pm is_j^y)$. This constitutes an

<figure><img src="image_4.png" /><figcaption>FIG. 3: (Color online) Main panel: Behaviour of $z \cdot C(\mathbf{r}, z)$ vs $r^2/z$, on a semi-log scale for $z = 0$ (×), 2 (●), 4 (■), 6 (▲), 8 (+), 10 (▼). Dashed line shows behaviour expected from Eq. (5). The system size is $L_\perp = 14$ and $L_z = 196$. Inset: Same data plotted as $C(\mathbf{r}, z)$ vs $r/L_\perp$ on a linear scale, showing that correlations vanish within (001) planes.</figcaption></figure>

enhancement of the symmetry at the critical point from U(1) to SU(2).

The quantum description can be employed to calculate correlation functions within a given sector in spin ice, by approximating the imaginary time direction as continuous and treating $\mathcal{H}_Q$ using harmonic spinwave theory. For a system of shape sufficiently anisotropic ($L_z \gg L_\perp^2$) to justify a ground state treatment of the 2+1 dimensional problem, and taking $z$ and $\mathbf{r}$ to denote distances along the [001] direction and in the (001) plane respectively, with $\rho$ a microscopic length, two-point correlations at $T_c$ are

$$ C(\mathbf{r}, z) \propto \frac{1}{z} \exp \left( - \frac{|\mathbf{r}|^2}{\rho z} \right). \qquad (5) $$

As the strings do not interact, this form reflects the string auto-correlations present in a random walk in two dimensions with propagation time $z$. It agrees well with results from simulations, as we show in Fig. 3.

A further consequence of the absence of interactions between strings is that the surface tension between oppositely magnetised domains vanishes and the domain wall width diverges, as $T_c$ is approached from below. We investigate this phenomenon in simulations by using weak, position-dependent magnetic fields to induce two domains with the interfaces between them lying on average in the x-z plane. Results are shown in Fig. 4. At $T_c^-$ the magnetisation profile $M(y)$ is a function only of $y/L_\perp$. Below $T_c$ a mean field treatment of $\mathcal{H}_Q$ gives a domain wall width $\ell_w^{-1} \propto \sqrt{t} \equiv (1 - T/T_c)^{1/2}$: this is confirmed in Fig. 4, up to a constant vanishing in the thermodynamic limit. It is striking to find broad domain walls in an Ising magnet; they may be detectable using small angle neutron scattering.

The natural way to generate the required degeneracy lifting in a magnetic compound is to apply uniaxial pressure along the [001] axis of a single crystal. We require
<figure><img src="image_5.png" /><figcaption>FIG. 4: Left: Magnetisation profile vs $y/L_⊥$ at $T = 0.999 T_c$ for $L_⊥ = \{6, 8, 10, 12, 16, 20\}$. Right: Domain wall width $ℓ_w$ vs $\sqrt{t}$ for $L_⊥ = 20$. The straight line is a guide to the eye. All samples have $L_z = L_⊥^2$.</figcaption></figure>

$\delta$ positive. Within the nearest neighbour approximation the effective coupling has two contributions, one from the long-range dipolar interactions [11, 12, 25], the other due to superexchange. While a uniaxial compression along [001] increases the former, the change of the latter is less clear as it depends on the evolution of orbital overlaps.

<figure><img src="image_6.png" /><figcaption>FIG. 5: Percentage change in magnetisation measured under uniaxial pressure of 10.5 kbar, $M(P, B)$, compared to that at atmospheric pressure, $M(0, B)$, at $T = 1.7$K as a function of applied field $B$ along the [001] axis (extracted from [26]).</figcaption></figure>

Results from the one experiment so far performed on Dy$_2$Ti$_2$O$_7$ under pressure [26] suggest that the effective coupling is indeed modified so as to produce $\delta > 0$. Pressure causes a 4% increase in $M$ at small field, but a 1% decrease at high field (Fig. 5). The increase in the zero field susceptibility is consistent with the lifting of the degeneracy in favour of the states with magnetization along the field axis, as required for the transition studied here (Fig. 1). The reduction of the saturated moment is expected as the crystal is squeezed along [001] because of a resulting tilt of the easy axes away from [001].

In fitting the parameter $\delta$, we account for demagnetisation effects by modelling the sample as a prolate ellipsoid with major axes given by the dimensions of the approximately rectangular parallelepiped specimen. For the interaction parameter of Dy$_2$Ti$_2$O$_7$ [27], the best fit, shown in Fig. 5, yields $T_c \sim 200$ mK.

Whereas this temperature is comfortably within reach of cryogenics, the dynamics of spin ice slows down greatly below $T_f \sim 600$ mK, as activated ice-rule violating defects disappear. Their presence thus appears indispensable in practice, even though they will lead to more standard second-order phase transition in the three dimensional Ising universality class [16, 19]. The width of its critical region will, however, vanish with defect concentration. Experiments should thus aim at the lowest dynamically accessible transition temperatures, with the target $T_c \sim 600$ mK only a factor of 3 away from existing experiments [26]. We hope our work will stimulate further experimental efforts to realise this unusual transition.

We thank K. Matsuhira for sharing results from Ref. [26], and D. Simon and S.T. Bramwell for useful discussions. We acknowledge financial support from ESF under Grant PESC/RNP/HFM (PCWH and LJ), from ANR under Grant 05-BLAN-0105 (LJ), and from EPSRC (JTC) under Grant EP/D050952/1.

[1] C. L. Henley, arXiv:0912.4531 (2009).

[2] F. Alet et al., Phys. Rev. Lett. 97, 030403 (2006).

[3] T. S. Pickles, T. E. Saunders, and J. T. Chalker, Europhys. Lett. 84, 36002 (2008).

[4] G. Chen et al., Phys. Rev. B 80, 045112 (2009).

[5] I. Mirebeau et al., Nature 420, 54 (2002).

[6] I. Mirebeau et al., Phys. Rev. B 74, 174414 (2006).

[7] M. E. Lines and A. M. Glass, *Principles and Applications of Ferroelectrics and Related Materials* (Clarendon, Oxford, 1977).

[8] J. C. Slater, J. Chem. Phys. 9, 16 (1941).

[9] E. H. Lieb, Phys. Rev. Lett. 19, 108 (1967).

[10] M. J. Harris et al., Phys. Rev. Lett. 79, 2554 (1997).

[11] M. J. P. Gingras and B. C. den Hertog, Can. J. Phys. 79, 1339 (2001).

[12] S. V. Isakov, R. Moessner and S. L. Sondhi, Phys. Rev. Lett. 95, 217201 (2005).

[13] S. V. Isakov et al., Phys. Rev. Lett. 93, 167204 (2004).

[14] C. L. Henley, Phys. Rev. B 71, 014424 (2005).

[15] S. V. Isakov et al., Phys. Rev. B 70, 104418 (2004).

[16] L. D. C. Jaubert, Ph.D. thesis, ENS de Lyon (2009).

[17] K. Binder, Rep. Prog. Phys. 50, 783 (1987).

[18] H. Takahashi, Proc. Phys. Math. Soc. Japan 23, 1096 (1941).

[19] Y. Takagi, J. Phys. Soc. Jpn 3, 273 (1948).

[20] J. F. Nagle, Comm. Math. Phys. 13, 62 (1969).

[21] J. D. Shore and D. J. Bukman, Phys. Rev. Lett. 72, 604 (1994).

[22] L. Benguigui, Phys. Rev. B 16, 1266 (1977).

[23] R. J. Baxter, *Exactly Solved Models in Statistical Mechanics* (Dover Publications, Mineola, 2007), Chap. 8.

[24] S. Powell and J. T. Chalker, Phys. Rev. B 78, 024422 (2008).

[25] R. Siddharthan et al., Phys. Rev. Lett. 83, 1854 (1999).

[26] M. Mito et al., J. Magn. Magn. Mater. 310, E432 (2007).

[27] B. C. den Hertog and M. J. P. Gingras, Phys. Rev. Lett. 84, 3430 (2000).