Thermodynamics and fractal dynamics of a nematic spin ice: A doubly frustrated pyrochlore Ising magnet

The Ising antiferromagnets on the triangular and on the pyrochlore lattices are two of the most iconic examples of magnetic frustration, paradigmatically illustrating many exotic properties such as emergent gauge fields, fractionalization, and topological order. In this paper, we show that the two instances of frustration can, remarkably, be combined in a single system, where they coexist without inducing conventional long-range ordering. Our results indicate that the system undergoes a first-order phase transition upon lowering the temperature, into a yet different frustrated phase that we characterize to exhibit nematic order. We argue that an extensive degeneracy survives down to zero temperature, at odds with a customary Pauling estimate. Dynamically, we find evidence of anomalous noise in the power spectral density, arising from an effectively fractal anisotropic motion of monopoles at low temperature.

Figure 1

Spin ice consists of a pyrochlore lattice with classical Ising spins that are constrained to point directly in or out of a tetrahedron $[{\sigma }_i=\pm 1$ in Eq. (1), with opposite convention on each of the two tetrahedral sublattices]. The ferromagnetic nearest-neighbor interactions generate an extensive number of ground states with two spins pointing in and two spins pointing out of every tetrahedron. These “2-in 2-out” requirements are known as the ice rules. In our model, Eq. (1), antiferromagnetic interactions between third-neighbor spins are also included in the direction perpendicular to the spin axis (conventionally referred to as $J_3^{\prime }$ coupling). Examples of these interactions are illustrated by the thick green lines, which also illustrate how these interactions form triangular planes cutting through the pyrochlore lattice. To avoid cluttering, these interactions are only displayed for one of the four spin sublattices—the others lie in the three symmetry-equivalent {111} planes. The black box indicates the size of a system with $L=1$ (i.e., the conventional 16-spin cubic unit cell of the pyrochlore lattice).

Figure 2

Illustration of the ${\beta }_{12}$ (orange) and and ${\beta }_{34}$ (cyan) spin chains, pointing along the [110] and [–110] directions, respectively. The numbers in black indicate the four spin sublattices. The black box indicates the size of a system with $L=1$ (i.e., the conventional 16-spin cubic unit cell of the pyrochlore lattice). Every spin belongs to one of these $\beta$ chains, and any state where all spins along each chain are aligned satisfies the ice rules. In addition, if the direction of these chains of aligned spins is chosen correctly (see main text), the energy of the third-neighbor interactions will also be minimized, generating a ground state of $\mathcal{H}$. Two alternative pairs of $\beta$ chains can be defined in the same way, ${\beta }_{13}$ and ${\beta }_{24}$ or ${\beta }_{14}$ and ${\beta }_{23}$.

Figure 3

Illustrative phase diagram of the Hamiltonian in Eq. (1). Both the high-temperature (cooperative) paramagnetic phase, shown in blue, and the low-temperature spin nematic phase, shown in red, have a spin ice regime at low temperature. The phases are separated by a phase transition at $T_c/J_1\approx 2.1{(J_3^{\prime }/J_1)}^{0.8}$, with the critical temperature determined in Figs. 5 and 6.

Figure 4

Spin correlations along the three chain directions ${\beta }_{12}$ (circles), ${\beta }_{13}$ (squares), and ${\beta }_{14}$ (triangles) for $J_3^{\prime }/J_1=0.2$ at temperatures $T/J_1=0.2<T_c/J_1$ (green) and $T/J_1=0.9>T_c/J_1$ (pink). The three directions correspond to $\mathbf{r}=(m,m,0), \mathbf{r}=(m,0,m)$, and $\mathbf{r}=(0,m,m)$ respectively, with integer $m$. Below $T_c$ one chain direction (${\beta }_{12}$ in this case) is chosen along which spins tend to align head to tail, resulting in large spin correlation along this direction. The staggered behavior of the spin correlations in the other directions, with nonzero correlations for even $m$, is a result of the antiferromagnetic interactions between chains of the same type (see Appendix pp5 for details).

Figure 5

Behavior of the magnitude of $\langle \mathrm{\Gamma }\rangle$, defined in Eq. (3). Going from left to right, the value of $J_3^{\prime }$ increases by factors of 4 between $J_3^{\prime }/J_1=0.05$ (orange) and $J_3^{\prime }/J_1=51.2$ (pink). The blue line corresponds to $J_3^{\prime }=0$ (i.e., conventional nearest-neighbor spin ice, where $\mathrm{\Gamma }$ vanishes at all temperatures). The black line corresponds to infinite $J_3^{\prime }$. $\langle \mathrm{\Gamma }\rangle$ acts as an order parameter, remaining zero up to the point where the system enters the nematic spin ice phase below the $J_3^{\prime }$-dependent critical temperature $T_c$. We define pragmatically $T_c$ as the highest temperature where $\left|\langle \mathrm{\Gamma }\rangle \right|>{10}^{-2}$, and we plot it as a function of $J_3^{\prime }$ in the inset (notice the logarithmic scale on both axes). The dotted black line $T_c/J_1=2.1{(J_3^{\prime }/J_1)}^{0.8}$ is a fit to the data.

Figure 6

Specific heat per spin for $\mathcal{H}$ with different values of $J_3^{\prime }/J_1$ (shown in the legend), computed for different system sizes: $L=3$ (dotted lines), $L=6$ (dashed lines), $L=9$ (dash-dotted lines), and $L=12$ (solid lines). The phase transition manifests itself as a peak in the specific heat, and the peak position depends only weakly on system size. The inset shows the critical temperature $T_c$, defined by the position of the specific heat peak for $L=12$, for each of the $J_3^{\prime }$ values. The black dotted line shows $T_c/J_1=2.1{(J_3^{\prime }/J_1)}^{0.8}$, demonstrating the good agreement with the behavior inferred from the order parameter in Fig. 5.

Figure 7

Histograms of the system energy at temperatures around the phase transition, for $J_3^{\prime }/J_1=0.2$ (left), $J_3^{\prime }/J_1=0.8$ (centre), and $J_3^{\prime }/J_1=3.2$ (right). The data were collected from 1000 independent simulations with $L=9$. From bottom to top (the curves are shifted for visualization purposes), the temperature is varied from $0.506J_1$ to $0.538J_1$ for $J_3^{\prime }/J_1=0.2$, from $1.650J_1$ to $1.698J_1$ for $J_3^{\prime }/J_1=0.8$, and from $4.578J_1$ to $4.994J_1$ for $J_3^{\prime }/J_1=3.2$. The temperature ranges were chosen to encompass the specific heat peak in each case, and the sample temperatures are distributed evenly in $logT$ within each range.

Figure 8

Binder cumulants for the real (filled circles and solid lines) and imaginary (unfilled circles and dashed lines) parts of the order parameter $\mathrm{\Gamma }$. The results are shown, from top to bottom, for $J_3^{\prime }/J_1=0.2, J_3^{\prime }/J_1=0.8$, and $J_3^{\prime }/J_1=3.2$, for temperatures around the phase transition. The Binder cumulants were computed from 100 independent simulations, for each of the three system sizes indicated in the legends.

Figure 9

The entropy computed by integrating the specific heat for a system of size $L=9$. The patterned black lines indicate entropy estimates of different models: the infinite temperature entropy per spin, $k_{B}ln\left(2\right)$, is shown as a dashed line; the dash-dotted line indicates Wannier's estimate of the entropy of the triangular lattice Ising antiferromagnet [2], 0.323066; and the dotted line indicates the Pauling entropy of conventional spin ice [24], $k_{B}ln(3/2)/2$. For the system size shown here, $L=9$, we observe a finite residual entropy at the lowest temperature of about $k_B/20$ per spin, for any $J_3^{\prime }/J_1>0$. This value is significantly higher than the (subextensive) entropy of the fully aligned ground states of $\mathcal{H}$, which tends to zero as $L^{-1}$ and is indicated by the grey line in the figure. These results are in agreement with a phase transition into an extensive subset of the ice states at low temperature.

Figure 10

The number density of magnetic monopoles (top) and of triangular lattice excitations (bottom) vs temperature for a range of $J_3^{\prime }/J_1$ values in a system of size $L=9$. Solid lines indicate simulations where the system was cooled from high temperature, whereas dashed lines indicate simulations begun from a randomly chosen ground state and then heated up. Vertical dotted lines of the corresponding color indicate the values of $T_c$ obtained from the order parameter (equiv., from the specific heat data).

Figure 11

PSD measured along the $x$ direction (top) and the $z$ direction (bottom) at the same temperature $T/J_1=0.4$, spanning different values of $J_3^{\prime }$. The system is in the nematic spin ice phase for $J_3^{\prime }/J_1>0.13$. Note that for $J_3^{\prime }/J_1\ge 0.4$ there appears to be no further change in the power-law behavior of the noise with variations in $J_3^{\prime }$. That the noise is already anomalous for $J_3^{\prime }/J_1=0.1$ suggests that the anomalous behavior appears already above the phase transition when $J_3^{\prime }/T\sim 1$ and the $J_3^{\prime }$ interactions begin to influence monopole motion. The noise continuously becomes more anomalous across the phase transition; once nematic correlations fully onset, the anomalous behavior remains unchanged irrespective of how deep one takes the system into the nematic phase by tuning $J_3^{\prime }$.

Figure 12

(Top) Average number of accessible sites for a monopole within chemical distance $n$, computed from 2000 independent configurations. The round markers show the behavior in the nematic phase when the monopole either makes all moves allowed by the ice rules (orange) or only moves that do not increase the energy (green). The clusters for bSM spin ice [16] (no $J_3^{\prime }$ interaction) and the critical percolation cluster on the diamond lattice are shown for comparison (dark blue squares and light blue crosses, respectively), see main text for details. The agreement between green round markers and light blue crosses suggests that monopoles appear to move on a fractal percolation cluster up to $n\sim 20$. (Bottom) Examples of monopole clusters in nematic spin ice when the monopole is either moving without accounting for the $J_3^{\prime }$ interactions (top row) or when it is only allowed to make moves that do not increase the energy (bottom row).

Figure 13

Mean-squared displacement of individual monopoles for a range of values of $J_3^{\prime }$. The incipient plateau appearing at long times is a finite-size effect (this is a system with $L=15$). $\langle Y^2\left(t\right)\rangle$ behaves like $\langle X^2\left(t\right)\rangle$, whereas $\langle Z^2\left(t\right)\rangle$ clearly grows more slowly, as one would expect from the behavior of the PSD. In these simulations a single monopole pair is created in each realization, and only one of the two monopoles is allowed to move (see Appendix pp3).

Figure 14

(a) Section of spin ice in a long-range ordered ground state of $\mathcal{H}$ with spins aligned along the ${\beta }_{12}$ (orange) and ${\beta }_{34}$ (cyan) chains. The $J_3^{\prime }$ interactions are shown for spins on sublattice 1 (green lines) and sublattice 2 (purple lines). If one projects the chains onto the plane normal to the chain direction (shown in grey for the ${\beta }_{12}$ chains), one is left with a triangular pattern as shown in panel (b). (b) Combined triangular pattern from the projection of the ${\beta }_{12}$ chains. A red dot indicates that the spins along this chain point out of the plane and a blue dot indicates that they point into the plane—one can think of each chain as a net Ising spin. Each bond on the combined pattern represents $4L{J}_3^{\prime }$ interactions in spin ice—including both the ones between spins on sublattice 1 [green lines in (a)] and between spins on sublattice 2 [purple lines in (a)]. Directing the chains so that each triangle on the triangular pattern has two blue and one red dots or one blue and two red dots therefore minimises the energy of all the $J_3^{\prime }$ planes formed by spins on sublattice 1 and 2.

Figure 15

Spin correlations along the ${\beta }_{12}$ chain directions (left panel), and along ${\beta }_{13}$ (right panel), from Monte Carlo simulations of systems of size $L=9$ with $J_3^{\prime }/J_1=0.8$ at temperatures near $T_c$ ($T_c/J_1\approx 1.7$). The coordinate system has been chosen so that the selected chain pairs at low temperature are ${\beta }_{12}$ and ${\beta }_{34}$. Alignment of the vector spins along a chain corresponds to Ising variables that follow the pattern $+,-,+,-,...$; this results in positive values of ${(-1)}^m\langle \sigma \left(\mathbf{0}\right)\sigma \left(\mathbf{r}\right)\rangle$ for $\mathbf{r}$ along the ${\beta }_{12}$ or ${\beta }_{34}$ directions. The staggered behavior of the correlators along other chain directions, an example of which is given in the right panel, is caused by the effective antiferromagnetic interaction between neighboring ${\beta }_{12}$ or ${\beta }_{34}$ chains.

Figure 16

Correlations of the Potts variable $\gamma$, defined in Eq. (2), along various symmetry directions. All plots are for $T/J_1=0.4$ (the temperature at which we measured the magnetic noise in Sec. 4), which is greater than $T_c$ for $J_3^{\prime }/J_1<0.1$. For $J_3^{\prime }/J_1\ge 0.8$ the data points overlap perfectly. The states have been chosen so that the selected chain pairs at low temperature are ${\beta }_{12}$ and ${\beta }_{34}$, explaining why correlations are stronger along the [110] and [1–10] directions. The numerical results show no difference between the [100], [010], and [001] directions, which is why only the first one is shown. Similarly, the [101] and [011] directions are equivalent. This is for system size $L=9$. The smaller values when $m$ is a multiple of 6 for certain directions is again a symptom of the antiferromagnetic interactions between chains of the same type and the triangular pattern that they form.

Figure 17

The densities of magnetic monopoles (blue) and triangular excitations (red) vs temperature for two values of $J_3^{\prime }$. Cooling and heating runs are shown by solid and dashed lines, respectively. In the top panel, with $J_3^{\prime }=0.2J_1$, the system falls out of equilibrium when cooled below the equilibrium transition temperature, indicated by the dotted grey line. In the bottom panel, with $J_3^{\prime }=0.4J_1$, there is no evidence of low-temperature out-of-equilibrium behavior (although some hysteresis is visible in the triangular excitation density). The dotted horizontal lines show the $1/V$ limit for the two types of excitations – below these lines there is on average less than one excitation of the respective type in the entire finite system in our simulations. The figure is for $L=9$.

Figure 18

PSD computed from the magnetization measured along the three different Cartesian axes. They are all at $T/J_1=0.4$ and with $J_3^{\prime }/J_1=0.1$ (top), 0.2 (middle), and 0.8 (bottom). Cooling runs are shown with solid lines and heating runs with dashed lines. The dotted black lines show the result of a fit with Eq. (6), with the extracted exponent shown in the legend. In the nematic spin ice phase, monopole transport along one axis (here chosen without loss of generality to be the $z$ axis) is reduced, as further discussed in Sec. 5 in the main text.