Monte Carlo study on the detection of classical order by disorder in real antiferromagnetic Ising pyrochlores

TL;DR

This study uses Monte Carlo simulations to explore the detectability of thermal order by disorder in real antiferromagnetic Ising pyrochlores, considering various perturbations and experimental methods.

Contribution

It introduces a flexible simulation approach to distinguish entropic order from interaction-driven order in real magnetic materials.

Findings

Conditions identified for observing entropic order effects.

Neutron scattering proposed as a detection method.

Analysis of effects of dipolar interactions and field misalignment.

Abstract

We use Monte Carlo simulations to evaluate the feasibility of detecting thermal order by disorder in real antiferromagnetic Ising pyrochlores, frustrated by a magnetic field applied in the $[110]$ direction. Building on an ideal system with only nearest-neighbour exchange interactions and a perfectly oriented field, we consider the effects of dipolar interactions and field misalignment. Our approach is special in that it relies more in the possibility to switch on and off the entropic drive towards order than in the absence of (or immunity to) a particular perturbation. It can then be applied, in principle, to other uncontrolled interactions expected to be naturally present in real {magnetic} materials. We establish the conditions under which entropic effects can be discerned from an interaction drive towards order, show how to use neutron scattering as a means to unveil this mechanism, and discuss possible materials where to test these ideas.

Equations

$$\mathcal{H}=$$

$$+D{r}_{\text{nn}}^3\sum _{i>j}\left[\frac{{\mathbf{S}}_i\cdot {\mathbf{S}}_j}{|{\mathbf{r}}_{ij}{|}^3}-\frac{3({\mathbf{S}}_i\cdot {\mathbf{r}}_{ij})({\mathbf{S}}_j\cdot {\mathbf{r}}_{ij})}{|{\mathbf{r}}_{ij}{|}^5}\right],$$

$$\frac{\mathcal{H}}{|J_{\text{eff}}|}=-\sum _{⟨ij⟩}{S}_i{S}_j-\frac{\sqrt{2}}{\sqrt{3}}\frac{h}{|J_{\text{eff}}|}\sum _{i\in \alpha }{\alpha }_i{S}_i+\frac{D}{|J_{\text{eff}}|}{\mathcal{H}}_{\text{dip}}^{r>r_{\text{nn}}},$$

$$\rho_{S}^{s}=\bigg{\langle}\frac{1}{N_{\textit{up}}}\left|\sum_{\nu}a_{\nu}^{s}\right|\bigg{\rangle}\,,$$

$$\rho_{S}^{d}=\bigg{\langle}\frac{1}{N_{\textit{up}}}\left|\sum_{\nu}a_{\nu}^{d}\right|\bigg{\rangle}\,,$$

$$I_{\text{SS}}(\mathbf{k})=\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^N⟨S_i{S}_j⟩\left({\hat{\mathbf{s}}}_i^{\perp }\cdot {\hat{\mathbf{s}}}_j^{\perp }\right)e^{i\mathbf{k}\cdot {\mathbf{r}}_{ij}},$$

$${\hat{\mathbf{s}}}_i^{\perp }={\hat{\mathbf{s}}}_i-\left({\hat{\mathbf{s}}}_i\cdot \frac{\mathbf{k}}{|\mathbf{k}|}\right)\frac{\mathbf{k}}{|\mathbf{k}|}.$$

$$I_{\text{QQ}}(\mathbf{k})=\frac{2}{N}\sum _{\eta =1}^{N/2}\sum _{\nu =1}^{N/2}⟨Q_{\eta }{Q}_{\nu }⟩e^{i\mathbf{k}\cdot {\mathbf{r}}_{\eta \nu }},$$

$$M=-\frac{{\mu }_0}{4\pi }\frac{\alpha Q^2}{r_d},$$

$$E^{\text{self}}=\left[J_{\text{eff}}-\frac{5}{3}D+\frac{4}{5}\left(1+\sqrt{\frac{2}{3}}\right)\frac{5}{3}D\right]{\left(\frac{r_d}{\mu }Q\right)}^2.$$

$$E_d=E_s-\frac{4}{\sqrt{6}}{h}_{\text{inf}},$$

Full text

I Introduction

In magnetism, the geometry of the spin lattice can be just an unimportant detail (as when considering critical fluctuations in a ferromagnet Cardy (1996)) or a crucial one, deciding whether magnetic order actually takes place or not. The latter case applies to geometrically-frustrated magnetic systems, characterised by the absence of spin order at low temperatures even without any quenched structural disorder Ramirez (1994); Moessner and Ramirez (2006). Prominent examples among these disordered systems are the magnetic equivalents of a liquid (spin liquids, correlated but disordered) Balents (2010); Lacroix et al. (2011) and of ice (spin ices) Bramwell and Gingras (2001); Diep (2013). They are home to very exotic quasiparticles Balents (2010); Fennell et al. (2014) and to the magnetic analogue of electric charges or “monopoles”, respectively. The monopoles interact through a Coulomb law and respond to a magnetic field like charges to an electric one Castelnovo et al. (2008). In further analogy to systems of electric charges, different states of monopole matter have been suggested or detected: monopole fluids Slobinsky et al. (2018); Borzi et al. (2013); Guruciaga et al. (2014), ionic-like crystals with staggered magnetic-charge density Castelnovo et al. (2008); Borzi et al. (2013); Brooks-Bartlett et al. (2014); Guruciaga et al. (2014); Jaubert (2015), or much more complicated structures Sazonov et al. (2012); Jaubert and Moessner (2015); Borzi et al. (2016); Udagawa et al. (2016); Rau and Gingras (2016). Much of these physics can be extended into two-dimensional artificial structures Xie et al. (2015); Farhan et al. (2019); Reichhardt et al. (2012), opening also to new possibilities.

Since strong spin correlations are operative, in principle even very small perturbations can break the massive degeneracy of some of the fluid-like (or even crystalline Borzi et al. (2013); Brooks-Bartlett et al. (2014)) ground states, forcing the appearance of long-range order at low-enough temperature. Adding to this, the fact that the ground state multiplicity in these frustrated systems is accidental opens up the possibility to other –more subtle– mechanisms of symmetry breaking. Seemingly against thermodynamics, even thermal disorder can be the driving force helping to develop an order parameter. Mirroring this apparent contradiction, the phenomenon was baptised as classical or thermal “order by disorder” Chalker (2011); Moessner and Chalker (1998); Shender and Holdsworth (1996); Bramwell et al. (1994) (OBD). It was first shown by Villain and collaborators in the so called generalised domino model Villain et al. (1980), and it occurs when there is a huge disproportion in the density of low-energy excitations associated with a particular ground state.

Theoretical models where disorder (including even structural defects, see for example Ref. Maryasin and Zhitomirsky, 2014) helps to develop a magnetic order parameter are quite interesting by themselves and have been studied for many years Villain et al. (1980); Henley (1989). However, a crucial question to ask when dealing with real systems is how to determine that the main cause of order is related to the entropy term ($-TS$) of the free energy $F=U_{0}-TS$, and not to an uncontrolled energy perturbation affecting the internal energy $U_{0}$ (see discussion on Ref. Savary et al., 2012). The strongest case in favour of actual OBD at the moment seems to be that of the XY pyrochlore $\mathrm{Er_{2}Ti_{2}O_{7}}$, where the classical ground state degeneracy is protected by symmetry Savary et al. (2012). The actual mechanism is still under discussion Rau et al. (2016), and quantum Savary et al. (2012); Zhitomirsky et al. (2012); Ross et al. (2014); Champion et al. (2003), classical McClarty et al. (2014) or composed Oitmaa et al. (2013) fluctuations could be the selection mechanism behind the observed order.

In this paper we explore in detail this crucial point of the detection of OBD in real materials. Following Ref. Guruciaga et al., 2016, we will see that antiferromagnetic Ising pyrochlores Sleight et al. (1974); Lhotel et al. (2015); Xu et al. (2015); Anand et al. (2015); Bertin et al. (2015) are a promising experimental scenario to detect OBD. While it is never possible to switch off every term affecting the unperturbed internal energy $U_{0}$ in a real crystal, the physics on these materials naturally suggest the reverse approach: the thermally-driven tendency towards magnetic charge order in the $-TS$ term can be detuned by means of an easily-controllable parameter. In this way, we manipulate the magnetic field to change the type of excitations occurring at fixed temperature, gauging the importance of the OBD mechanism by its effect on the staggered charge density. Using Monte Carlo simulations, we apply this method to the particular case of dipolar interactions. These, aside from being expected in real materials, allow understanding the perturbation as an attraction among Coulomb charges and thus favouring charge order. In principle, the same ideas could be applied to any other perturbation terms affecting $U_{0}$ (e.g., second or further neighbours exchange interactions), or to other systems presenting field-tuned OBD (e.g., artificial spin ices Guruciaga et al. (2016), where dipolar interactions are unavoidable). We also study the effect of a misaligned magnetic field $\mathbf{B}$. This case is special since the control parameter used to put in evidence the OBD mechanism is (in addition to this) a force driving the order; in other words, B changes $S(T,\mathbf{B})$ in a nontrivial way but also perturbs $U_{0}$. With the intention to bridge the gap between theory and experiment, we suggest materials where OBD could be detected using these ideas and how to measure the effect, concentrating in the particular case of neutron scattering techniques.

I.1 System and model

The pyrochlore lattice consists of corner-sharing tetrahedra, the centres of which form a diamond lattice (Fig. 1(a)). Classical Ising magnetic moments ${\bm{\mu}}_{i}=\mu{\mathbf{S}}_{i}=\mu S_{i}\hat{\mathbf{s}}_{i}$ sit on the vertices of the tetrahedra, with the quantisation directions $\hat{\mathbf{s}}_{i}$ along the local $\langle 111\rangle$ directions. The pseudospins $S_{i}=\pm 1$ indicate if the magnetic moments point outwards ($+1$) or inwards ($-1$) of “up” tetrahedra (coloured ones in Fig. 1(a)). The Zeeman coupling with a magnetic field $\mathbf{B}$, as well as the magnetic interactions of exchange and dipolar origin, of strengths $J$ and $D$, respectively, are accounted for by the Hamiltonian

[TABLE]

where $\langle ij\rangle$ means the sum is carried over nearest neighbours, $r_{\textit{nn}}$ is the pyrochlore nearest-neighbour distance and $D=\mu_{0}\mu^{2}/(4\pi r_{\textit{nn}}^{3})$. It can be easily seen in Fig. 1 that the $N$ spins may be separated into two groups according to how they couple to a field $\mathbf{B}\parallel[110]$: while $\beta$ spins (blue) are orthogonal to it, $\alpha$ spins (yellow) have a nonzero projection ${\hat{\mathbf{s}}}_{i}\cdot{\mathbf{B}}=\alpha_{i}\sqrt{2/3}\,B$ with $\alpha_{i}=\pm 1$ Hiroi et al. (2003). Using these definitions, and referring all energy contributions to the nearest-neighbour interactions term, the Hamiltonian can be rewritten as

[TABLE]

with $J_{\textit{eff}}=J/3+5D/3$ (assumed $<0$ throughout this text), and $h=\mu B$, with the second sum running over $\alpha$ spins only. $\mathcal{H}_{\textit{dip}}^{r>r_{\textit{nn}}}$ is unitless, and encompasses all dipolar contributions beyond nearest neighbours. The ferromagnetic version of Eq. (2) corresponds to the well-known dipolar spin-ice model Melko and Gingras (2004). In our work we will focus on a case in which the antiferromagnetic nearest-neighbour term ($J_{\textit{eff}}<0$) modified by the Zeeman energy dominates. As shown analytically and numerically in Ref. 37, Eq. (2) with $D/|J_{\textit{eff}}|=0$ conducts to physics similar to Villain’s domino model Villain et al. (1980), but its field-tuned OBD may allow for an experimental contrast of the phenomenon. Different from this recent work, the nearest neighbor Hamiltonian is now perturbed by dipolar interactions $\mathcal{H}_{\textit{dip}}^{r>r_{\textit{nn}}}$, with a relative intensity controlled by the (assumed small) parameter $D/|J_{\textit{eff}}|$.

Within the monopole picture Castelnovo et al. (2008), the diverse spin configurations can be mapped to arrays of magnetic charges that live in the dual diamond lattice, with a lattice constant $r_{d}=\sqrt{3/2}\,r_{\textit{nn}}$. While tetrahedra with two spins pointing in and two out are considered neutral, those with three spins in and one out, or viceversa, have a positive or negative single charge $Q_{s}=2\mu/r_{d}$ sitting in their centres (small spheres in Figs. 1(a)-(c)). Moreover, positive or negative double charges (big spheres in Figs. 1(a) and 1(b)), with $Q_{d}=2Q_{s}$, can be found in tetrahedra with all spins pointing in or all spins pointing out, respectively. It will be useful to define the density of single monopoles $\rho^{s}$, which has a value of $1$ when there is exactly one single monopole per tetrahedron; analogously, one can define the density of double monopoles $\rho^{d}$ and neutral sites $\rho^{n}$. 111Note that, differently from the staggered charge density defined below, these are number densities, irrespective of the monopole charges (which always average to zero). Having transformed a configuration of spins into one of monopoles, dipolar interactions can in turn be approximately described as Coulombian forces between monopoles, as was shown by Castelnovo and collaborators Castelnovo et al. (2008). Since these forces promote the development of ordered phases of charges of alternating sign with the structure of zinc blende Borzi et al. (2013); Brooks-Bartlett et al. (2014); Guruciaga et al. (2014), we define here two order parameters: the single-monopole (double-monopole) staggered charge density $\rho_{S}^{s}$ ($\rho_{S}^{d}$) as the thermal average of the modulus of the total magnetic charge due to single (double) charges in up tetrahedra, normalised so that full order corresponds to a value of $1$. Mathematically, the single-monopole staggered density can be written

[TABLE]

where the Greek index $\nu$ denotes a sum carried over the diamond-lattice sites corresponding to the centres of up tetrahedra (the total number of which is $N_{\textit{up}}=N/4$), and $a_{\nu}^{s}$ is equal to $1$ ($-1$) if there is a positive (negative) single monopole in site $\nu$, and equal to [math] otherwise. Analogously, the double monopole staggered density is defined as

[TABLE]

with $a_{\nu}^{d}$ equal to $1$ ($-1$) if a positive (negative) double monopole sits in site $\nu$ of the up tetrahedra sublattice, and equal to [math] otherwise. According to these definitions, having for example $\rho_{S}^{s}<1$ can imply both a lattice full of partially disordered single charges, as well as one with a low density of single charges but with a perfect staggered order. The measurement of the density of single and double charges erases any ambiguity.

I.2 Simulation details

We performed Monte Carlo simulations with the Metropolis algorithm and single-spin-flip dynamics. In order to implement Eq. (2) in the algorithm, we used Ewald summations to take into account long-range interactions. Different materials were explored by changing the adimensional parameters in the Hamiltonian. In practice, we varied the value of $D$ (for example, increasing it to consider a material with a big $\mu$) modifying $J$ accordingly so as to keep $J_{\textit{eff}}$ constant. $B$ was always varied freely, from which we obtained the product $h=\mu B$.

The conventional unit cell of the pyrochlore lattice (Fig. 1(a)) consists of $16$ spins, and we simulated cubic systems of $L^{3}$ cells with periodic boundary conditions. As an example, a system with $L=6$ and dipolar interactions took $5\times 10^{4}$ Monte Carlo steps to reach the equilibrium, and then $2.5\times 10^{4}$ steps were needed to calculate the averages at each value of temperature and applied magnetic field. In turn, the results were averaged over $30$ independent runs.

We used spin configurations obtained from Monte Carlo runs to simulate the spin structure factor given by

[TABLE]

where $N$ is the number of sites in the pyrochlore lattice, $\langle S_{i}S_{j}\rangle$ is the thermal average of the correlation between pseudospins at sites $i$ and $j$, and $\hat{\mathbf{s}}_{i}^{\perp}$ is the component of the quantisation direction at site $i$ perpendicular to the scattering wave vector $\bm{k}$, calculated as

[TABLE]

We also simulated the charge-charge correlation function by following the expression

[TABLE]

where $N/2$ is the number of tetrahedra, $Q_{\eta}$ represents the topological charge at site $\eta$ of the diamond lattice, and ${\bm{r}}_{\eta\nu}$ is the distance between monopoles. Both the spin structure factor and the charge-charge correlation function were obtained by averaging over sets composed of $150-200$ configurations for $L=8$.

II Effect of dipolar interactions

II.1 Ground state and low-energy excitations

While the phenomenon we want to investigate necessarily happens at finite temperature, we first need to study the system at the lowest temperatures in order to determine the magnetic-field region where it would be observable. The unperturbed situation, where Eq. (2) is limited to nearest-neighbour interactions ($D/|J_{\textit{eff}}|=0$), has been studied before in Ref. Guruciaga et al., 2016. In the absence of magnetic field, the ground state is a crystal of double charges alternating in sign, making a zinc blende structure with $\rho_{S}^{d}=1$ Guruciaga et al. (2014). Having no magnetic moment, double monopoles are unstable under any sufficiently strong magnetic field. Fig. 2(b) shows this for the relevant situation were the field is along the $[110]$ direction. We computed the change in energy $e(h)$ of a single tetrahedron in the nearest-neighbour approximation due to the increasing field $h$ in the Zeeman term. We classified spin configurations as neutral ($n$) sites, and double ($d$) and single ($s$) monopoles. Above an inferior field $h_{\textit{inf}}$ the $\alpha$-spin chains result completely polarised, and two particular single monopoles (one of each sign) are the stable configurations in each tetrahedron. Construction constraints involving the underlying spins impose antiferromagnetic order within $\beta$-spin chains, and they in turn form charge-ordered chains running perpendicular to the field (see Fig. 1(c)). The polarised $\alpha$ spins decouple each of these ordered $\beta$ chains from the rest. The ground state turns then into a $\rho^{s}=1$ disordered array with subextensive residual entropy Guruciaga et al. (2016), characterised by $\rho_{S}^{d}=0$ (there are no double charges) and $\rho_{S}^{s}=0$ (with $D=0$ there is no driving force at $T=0$ to impose coherence between the charge-alternating single-monopole $\beta$ chains).

Although there are no further changes in the ground state on increasing field, it is important for the study of OBD to establish the field range where double monopoles are the lowest-lying energy excitations out of the single-monopole ground state. It has been shown before that, due to construction constraints involving the spins that make the magnetic charges, these double monopoles favour pair monopole correlations even in the absence of dipolar interactions Guruciaga et al. (2014); Slobinsky et al. (2018), while neutral sites disfavour them. Applied to the present situation, this general argument has important consequences. Imagine neutral sites are so energetically costly compared with double charges that their density is negligible at low temperatures. It is then easy to check that the only way to promote low energy excitations —i.e., double monopoles— is by turning an $\alpha$ spin linking plus and minus single charges in consecutive charge-ordered chains (see Fig. 1(b) and (c)). As shown numerically and analytically in Ref. Guruciaga et al., 2016, the proliferation of these double charges at finite temperature drives the appearance of long-range charge order by the coherent assembly of adjacent $\beta$ chains Guruciaga et al. (2016). Fig. 2(b) allows us to extract the field range $[h_{\textit{inf}},h_{\textit{sup}}]$ (black dotted lines) where OBD was shown to be possible for the particular case $D=0$.

The inclusion of dipolar interactions ($D\neq 0$) is made conceptually and even quantitatively simple within the monopole picture. We expect the zinc-blende crystal of double charges to become more stable (through an extra term proportional to $-4DQ_{s}^{\,2}$) and thus to extend to higher fields. In turn, the disordered ground state of single monopoles should now be ordered, stabilised by a term proportional to $-DQ_{s}^{\,2}$. This is why —if dipolar energy was the main ordering force— we would expect the charge-ordered state of single monopoles to be destroyed at a critical temperature scaling with $D$ (i.e., comparatively small within the context of this study). This way of reasoning is what will allow us to observe indirectly OBD effects even in the presence of such perturbations: we will choose to work at a “high” temperature $T\approx|J_{\textit{eff}}|/k_{B}$ such that the effect of perturbations is small, with OBD as the dominating driving force towards order. As we will see in Sec. II.2, the need to work at a relatively high temperature will make the upper limit $h_{\textit{sup}}$ much less significant than $h_{\textit{inf}}$.

In yet another example of its simplicity and beauty, we estimate now the field $h_{\textit{inf}}$ using the monopole picture. The Madelung energy per pair of particles of a crystal of magnetic charges $Q$ and $-Q$ in the zinc-blende structure is given by

[TABLE]

with $\alpha=1.6381$ Ashcroft and Mermin (1976) and $Q$ a function of $D$ through its link to the magnetic moment $\mu=\sqrt{({4\pi}/{\mu_{0}})r_{\textit{nn}}^{3}\,D}$. Their self-interaction energy Castelnovo et al. (2008), in turn, is

[TABLE]

The factor on the right is a constant taking the value $4$ or $8$ (depending on whether the crystal is made of single or double charges, respectively), rendering $E^{{\textit{self}}}$ a function of $D$ and $J_{\textit{eff}}$ only. Taking this into account, $h_{\textit{inf}}$ can be calculated by simply considering the Zeeman energy of the configurations favoured by the field and solving

[TABLE]

where the energy $E_{s}=M_{s}+E_{s}^{{\textit{self}}}$ related to the single-monopole crystal is equal to the sum of Eqs. (8) and (9) evaluated at $Q_{s}$. Similarly, $E_{d}=M_{d}+E_{d}^{{\textit{self}}}$ is the energy of the double monopole crystal. We show the analytical predictions of the monopole picture for different values of $D$ in Fig. 2(a). We include also the results for $h_{\textit{inf}}$ and $h_{\textit{sup}}$ obtained by simulating small spin systems ($L=2$) at very low temperatures (down to $k_{B}T/|J_{\textit{eff}}|=0.18$) 222Although the lowest-possible temperatures are desirable in order to determine the ground-state transition field $h_{\textit{inf}}$, compromises have to be made not only to ensure equilibrium but also to resolve $h_{\textit{sup}}$, which involves excitations and is thus extremely subtle. and extrapolating to $T=0$ the fields where the number densities of the different tetrahedron configurations of interest cross. More explicitly, we looked for the condition $\rho^{s}=\rho^{d}$ to estimate $h_{\textit{inf}}$ (see Sec. II.2.1), and $\rho^{n}=\rho^{d}$ for $h_{\textit{sup}}$.333We will see in the next section that the density of double monopoles could be affected in a non-trivial way by very peculiar finite-size effects that induce disorder at very low temperatures for $D=0$ Guruciaga et al. (2016). However, the system remains ordered even at $T=0$ for $D\neq 0$, which is the important case here. Fig. 2(a) shows that, in spite of the approximations involved, both estimations for $h_{\textit{inf}}$ give very similar results.

II.2 Finite temperature

II.2.1 Nearest-neighbours interactions ($D=0$)

The monopole picture implies that there is no energetic drive towards charge order for $D=0$. In spite of this, at low fields and temperatures we expect to find a crystal of double charges with the zinc-blende structure. This “all-spins-in–all-spins-out” structure is stabilised not by charge-charge interactions, but due to the constraints imposed by the spins that make these charges Guruciaga et al. (2014). For $h\gtrsim h_{\textit{inf}}$ (which stabilises a disordered ground state of single monopoles) the situation differs from what it would be trivially expected. As first shown in Ref. Guruciaga et al., 2016, charge disorder should occur for $T$ strictly zero as predicted but, as we explained before and will see below, the proliferation of double charges favours an array of staggered magnetic charges at finite temperature. This array transitions to a disordered state only for $k_{B}T/|J_{\textit{eff}}|>1$.

Figs. 3, 4, and 5 show the results of Monte Carlo simulations for a field parallel to $[110]$ and $D=0$ (black symbols and lines), in a system with $L=6$. We can see in Fig. 3 that the number density of the different monopoles in a tetrahedron evolves with increasing magnetic field at $k_{B}T/|J_{\textit{eff}}|=0.9$ from a majority of double monopoles to one of single charges, as expected. At the highest fields, neutral sites are the main excitation. Due to the energy scales involved we expect to have at these temperatures only: i- positive and negative double monopoles, ii- the two types of single monopoles with magnetic moment direction favoured by the magnetic field, and iii- the same for neutral sites (see continuous lines in Fig. 2(b)). This considerably helps the analysis; for example, to understand that at the transition field $h_{\textit{inf}}$ the density of double and single charges is $\rho^{d}\approx\rho^{s}\approx 0.5$ (curves crossing in Fig. 3). In the same terms, it may be expected that the curves for the double-monopole and neutral-sites densities cross very near $h_{\textit{sup}}$. However, against intuition, we can see them in this last case crossing at much smaller fields than that. While we will be able to explain this when analysing the order parameter behaviour, we can advance now that this rapid proliferation of neutral sites at low fields implies the destruction of OBD at $h$ much smaller than $h_{\textit{sup}}$.

The relevant order parameters (the staggered charge density of double and single monopoles) are plotted in Figs. 4(a) and 4(b) as functions of field at the same working temperature $k_{B}T/|J_{\textit{eff}}|=0.9$. The relatively big value of $\rho_{S}^{d}$ in Fig. 4(b) at low fields reflects the zinc-blende structure of double monopoles, gradually fading with increasing $h$. The fact that $\rho_{S}^{d}\neq 0$ for $h\gtrsim h_{\textit{inf}}$ reveals that a considerable number of ordered double monopoles (created by thermal fluctuations) remain even above this threshold. The persistence of double-monopole order above $h_{\textit{inf}}$ and the continuous aspect of the curves on traversing the field threshold makes it tempting to think that these measurements —performed on increasing field— are affected by the past history, retaining at high fields some of the zinc-blende structure which characterises order at lower ones. This is not the case: the same results are obtained in field-decreasing runs at fixed $T$ (not shown), and at fixed field and decreasing temperature (Fig. 5(a)). As mentioned before and made clear in this same figure, although the system should be disordered at zero temperature, charge order exists at finite temperature below the scale of $|J_{\textit{eff}}|/k_{B}$.

The behaviour of $\rho_{S}^{s}$ seems to be complementary to that of $\rho_{S}^{d}$, with an increasing amount of single monopole order at fields above $h_{\textit{inf}}$. The establishment of staggered single-charge order at this finite $T$ above $h_{\textit{inf}}$ has been analytically and numerically studied as a textbook case of field-tuned OBD, triggered by a diverging peak in the density of states (see Fig. 3 in Ref. Guruciaga et al., 2016) associated with this particular monopole array. Regarding the single charge order at fields below $h_{\textit{inf}}$, it can be naturally interpreted in terms of the underlying “all-in–all-out” double-monopole crystal lattice, to which single charges are the lowest energy excitations (Fig. 2(b)). The converse is true for the observed finite value of $\rho_{S}^{d}$ for $h\gtrsim h_{\textit{inf}}$.

In relation to this, we owe an explanation for the reduced value of $\rho^{d}$ with respect to $\rho^{n}$ at moderate fields $h\lesssim h_{\textit{sup}}$, observed in Fig. 3. At very low temperatures neutral sites are almost banned by their high energy cost, while the excitation of the allowed quasiparticles (double monopoles) imply the existence of staggered single-charge order Guruciaga et al. (2016). At higher temperatures, the arousal of neutral charges help to decorrelate the $\beta$ chains, and the staggered order is lost (Fig. 5(a)). The absence of order among consecutive $\beta$ chains makes now less probable to excite double monopoles (which are mainly created by flipping $\alpha$ spins that link staggered-ordered $\beta$ chains (see Fig. 1(c) and Ref. Guruciaga et al., 2016)). Since neutral sites are not so constrained (i.e., they do not need the coherent order of consecutive $\beta$ chains to appear) this explains the crossing of their density curves at finite $T$ below $h_{\textit{sup}}$ observed in Fig. 3(a). To illustrate this phenomenon, Fig. 5(b) shows the density of monopoles, now at fixed $h$ and as a function of temperature (black lines and symbols). We would naively have expected a density of double monopoles monotonously increasing with temperature. We can see that while $\rho_{S}^{s}$ decreases near the transition at $k_{B}T\approx|J_{\textit{eff}}|$ (Fig. 5(a)) the number of double monopoles starts to decrease, producing a maximum in its density within the staggered-order region. This finite temperature phenomenon implies no contradiction with the method we used to find $h_{\textit{sup}}$, performed in the limit of zero temperature.

The abrupt fall on $\rho_{S}^{s}(h)$ at fixed temperature on increasing fields observed in Fig. 4(a) is one of the main results we have to show in this work. It relates to the switching-off of the entropic drive towards order (the only one present in the current situation) produced by the proliferation of neutral sites Guruciaga et al. (2014); Slobinsky et al. (2018). A key idea is that both the thermally-induced single-charge correlations produced by double monopoles at low fields, and the decorrelating effect of neutral sites that causes a drop in $\rho_{S}^{s}$ at higher ones, will still be operative if dipolar interactions or other perturbations are present. As will be shown in the next subsection, this fact can be used to mark the presence of OBD effects.

II.2.2 Dipolar interactions (finite $D$)

Figs. 4(a) and 4(b) present $\rho_{S}^{s}$ and $\rho_{S}^{d}$ as functions of the field when long-range dipolar interactions are included. These curves were measured in the same conditions as the nearest-neighbour case. It can be seen that as $D$ grows, an increasing field $h$ is needed to destabilise the staggered order of double monopoles in favour of an order of single monopoles, as predicted by Fig. 2(a). The shape of the curves $\rho_{S}^{s}(h)$ allows us to separate them (at this particular working temperature) in two distinct sets: values of $D$ smaller and bigger than $D/|J_{\textit{eff}}|=0.13$.

The group of smaller $D$ is the relevant one concerning the main objective of this paper. Within this set, the temperature $T=0.9|J_{\textit{eff}}|/k_{B}$ is so high that Coulomb forces among monopoles cannot hold the staggered charge array. The curves resemble then those obtained without dipolar long-range interactions, showing that the main effect of a finite $D$ is to shift $h_{{\textit{inf}}}$, but preserving the same physics at this temperature. The abrupt fall of the order parameter down to $\rho_{S}^{s}=0$ on increasing $h$, which persists in the thermodynamic limit (Fig. 6), mirrors the switching-off of the OBD mechanism and serves as an indication of its relevance in the stabilisation of order.

A non-trivial behaviour associated with OBD is exemplified by Fig. 5(a). There, the order parameter $\rho_{S}^{s}$ is compared for $D=0$ (black curves) and $D/|J_{\textit{eff}}|=0.06$ (blue), now as a function of temperature at $h/|J_{\textit{eff}}|=7.86$ and $8.16$, respectively, where the corresponding $\rho_{S}^{s}(h)$ in Fig. 4(a) has a maximum. Diverting from the blue curve, $\rho_{S}^{s}$ for $D=0$ tends to zero at very low temperatures (and, evidently, for a constant lattice size $L$). This has been shown in Ref. Guruciaga et al., 2016 to be a very peculiar finite-size effect characterising this system (on increasing the size of the lattice the ordering trend continues on decreasing temperature, with disorder prevailing only at $T=0$) that could be used as a method to identify OBD in artificial two-dimensional materials. On the other hand, $\rho_{S}^{s}$ continues growing when the system is cooled down at finite $D$: on reducing thermal fluctuations such that $k_{B}T/D\to 0$, dipolar interactions are able to promote monopole order. Regarding this usual measure of the relative value of thermal fluctuations compared with the tendency towards order (here, $k_{B}T/D$), there is a further aspect to note. Near the order-disorder transition ($k_{B}T/|J_{\textit{eff}}|\approx 1$ for both cases in Fig. 5(a)) this measure is quite big, $k_{B}T/D\approx 20$ for $D/|J_{\textit{eff}}|=0.06$ (of the order of $100$ for $D/|J_{\textit{eff}}|=0.01$), while it diverges for $D=0$. This makes us wonder once again about how much against intuition it is that charge order could be sustained in the ideal case in spite of the total absence of charge attraction and repulsion.

For the second group ($D/|J_{\textit{eff}}|\gtrsim 0.13$) the shape of the curves in Fig. 4(a) starts to distort considerably. Their most salient feature is that there is no drop in the order parameter $\rho_{S}^{s}$ at high fields (where the corresponding single-charge densities are essentially saturated, as observed for the orange curve in Fig. 3). This indicates that the Coulomb forces can now maintain the order of charges at this temperature. It can also be noted that $\rho^{s}\approx 1$ at high fields: the concentration of doubles is made unstable by the Zeeman energy, while neutral sites are disfavoured by the Coulomb term hidden in Eq. (2).

Being $\rho_{S}^{s}$ and $\rho_{S}^{d}$ quantities generally inaccessible to most laboratory-based experiments 444The measurement of magnetocapacitance may be an interesting exception to this statement, allowing an indirect approach to the formation of staggered monopole order (see Refs. Katsufuji and Takagi, 2004; Saito et al., 2005; Khomskii, 2012)., an interesting probe to search for classical OBD in these magnetic systems would be neutron scattering. Figs. 7(a)-(c) show the simulated spin structure factor (Eq. (5)) for a system with $D/|J_{\textit{eff}}|=0.06$ and $L=8$ at the working temperature $k_{B}T/|J_{\textit{eff}}|=0.9$. We have chosen three values of the field, representative of the three types of order (or lack thereof) in the system. The patterns reflect the changes in magnetic order as the field turns the “all-in–all-out” structure (Fig. 7(a)) into a polarised crystal of single monopoles (Fig. 7(b)), and then into a disordered array of ordered $\beta$-spin chains (Fig. 7(c)). The evolution of the $(220)$ peak can be used to give a quantitative measure of three-dimensional order (i.e. between different $\beta$ chains). Supporting this, Figs. 7(d)-(f) show the charge-charge correlation function defined in Eq. (7) (which cannot be directly measured using neutrons, but illustrates well the parallel situation in the monopole framework). We can see that the sharp Bragg peaks associated with a crystal of monopoles change their intensity when passing from double to single charge order, and finally vanish for $h\approx h_{\textit{sup}}$, where the proliferation of neutral excitations promotes the disordering of the system.

There exist many “all-in–all-out” materials in nature, some of which could be used for an experimental contrast against these ideas Sleight et al. (1974); Lhotel et al. (2015); Xu et al. (2015); Anand et al. (2015); Bertin et al. (2015). However, it is important to bear in mind that each compound has its singularities, differing from the ideal case we study in its own way. $\mathrm{Nd_{2}Hf_{2}O_{7}}$ and $\mathrm{Nd_{2}Zr_{2}O_{7}}$, for example, are pyrochlores with a big Ising anisotropy. They have relatively small Néel temperatures ($T_{N}\lesssim 0.5\text{ K}$) and big magnetic moments on the $\mathrm{Nd}$ ion ($\approx 2.5\,\mu_{B}$), reduced by quantum fluctuations Lhotel et al. (2015); Anand et al. (2015); Xu et al. (2015); Anand et al. (2017). While $\mathrm{Nd_{2}Zr_{2}O_{7}}$ is a good candidate to test the behaviour of our system on the high $D/|J_{\textit{eff}}|$ side Xu et al. (2015), $\mathrm{Nd_{2}Hf_{2}O_{7}}$ has $D/|J_{\textit{eff}}|\approx 0.11$ Anand et al. (2015), slightly into the region where a decrease in order with increasing field would be expected. Finding an “all-in–all-out” compound with very small $D/|J_{\textit{eff}}|$ ratio to conclusively detect thermal OBD is, to our knowledge, much harder. The compound $\mathrm{Cd_{2}Os_{2}O_{7}}$, with quite a big ordering temperature (above $220\text{ K}$) Yamaura et al. (2012), could be a starting point, although it is clear that the physics implied (which has been studied since 1974 Sleight et al. (1974)) transcends that of a simple system with localised Ising moments.

III Effect of field misalignment

While it is pretty easy in a numerical simulation, applying a field perfectly oriented in a particular crystallographic axis is something practically impossible in the experimental realm. This fact strongly impacts on the possibility of detecting the OBD phenomenon in the way we propose here. A deviation of the field within the $x-y$ plane facilitates the formation of neutral tetrahedra without breaking the symmetry between the two FCC sublattices making the diamond structure. Therefore, it preserves spontaneous OBD but lowers the disordering temperature. On the other hand, a nonzero component of $\mathbf{B}$ in the $[001]$ direction explicitly breaks the relevant symmetry, and thus favours the occurrence of single-charge order. It is easy to understand this by noting that now the field has a larger component parallel to $[111]$, which stabilises the order of positive single charges in up tetrahedra and negative in down tetrahedra. 100100footnotetext: The $[110]$ direction has the special feature of being just midway between $[111]$ and $[11\overline{1}]$; the latter promotes the opposite order. Note (100); Sakakibara et al. (2003)

As mentioned in the introduction, this case is different from other perturbation terms, since the perturbation coupled with the order parameter grows in parallel with the magnetic field $h$, which is the instrument we use to detect the entropic drive towards order. To assess the significance of this effect at finite $T$, we consider the nearest-neighbour model (Eq. (I.1) with $D=0$) under a magnetic field $\mathbf{B}$ of strength $h/\mu$ tilted an angle $\theta$ from the $[110]$ direction towards $[001]$ (i.e., $\theta$ represents an anticlockwise rotation of that magnitude around $\hat{\mathbf{c}}\propto[1\overline{1}0]=[110]\times[001]$, see Figs. 1(a) and 1(b)). As shown in Fig. 8, the presence of $h_{001}=h\sin\theta\neq 0$ splits the low-lying single-charge energy levels associated to each type of tetrahedron (up or down), opening a gap proportional to $h_{001}$. For any value of $\theta\neq 0$, the single-charge ground state is strictly composed of one particular configuration and thus ordered. However, for small-enough values of $\theta$, the gap is effectively nonexistent in the range of $h$ and $T$ in which we are interested throughout this work, and the misaligned system can still show a signature of OBD.

Fig. 9 shows the relevant order parameter $\rho_{S}^{s}$ as a function of the modulus of the field, for different values of the tilting angle $\theta$. This figure is analogous to Fig. 4(a). It was measured at the same temperature and for the same system size, but with a different perturbation term. In black, we include as a reference the curve with a perfectly aligned field; given that $D=0$, it is the same black curve as the one in Fig. 4(a). The general shapes of the curves are similar to those encountered before. However, for all finite values of $\theta$, we observe at high $h$ a finite (and indeed, slightly increasing with field) value of $\rho_{S}^{s}$. As cautioned before, this marks an important difference with the previous perturbation: the system remains charge-ordered at high fields in any real situation.

Striking the main target we have fixed for this paper, the sudden drop observed in $\rho_{S}^{s}$ at very small angles (well below $1^{\circ}$) can only be explained by the weakening of the OBD effect that holds the order of single monopoles at intermediate fields even for $\theta=0$. Again, the observation of this drop is an indirect proof of the OBD mechanism being operative. The need of such an accuracy in the field angle control certainly complicates the observation of the effect, but it should not be considered as an experimental impediment: triple-axis vector magnets (as, for example, in Ref. Bruin et al., 2013), high-precision goniometers (like in Ref. Sazonov et al., 2010) or in-situ sample rotators Kenzelmann (2019) can be used to correct misalignments below the needed accuracy in a neutron scattering experiment.

IV Discussion

As advertised before, in our system the dipolar perturbation term that favours order remains constant when $h$ increases; this allows to recognise the effect in the order parameter of the field-tuned entropic contribution. However, it may be important to consider that bigger fields favour the introduction of neutral sites at the expense of double monopoles, not only weakening the thermal drive towards order, but also slightly diminishing the effectiveness of the Coulomb interactions. For small values of $D$ ($D/|J_{\textit{eff}}|<0.03$) the differences observed in monopole densities with respect to $D=0$ (not shown) are negligible, making evident that the OBD effect is the main driving force. For $D/|J_{\textit{eff}}|=0.06$, the field at which charge disorder occurs is somewhat bigger to that of $D=0$ (dashed lines in Fig. 3). Although the density of neutral sites at the transition fields is quite comparable, the density of double monopoles for $D/|J_{\textit{eff}}|=0.06$ is about half of the value it takes at the transition for $D=0$. This signals that, although dipolar interactions are not enough to hold the single-monopole crystal together at high fields, they are helping to sustain it near the order-disorder transition.

In a final note we will refer to the energy scales of the relevant perturbations needed to sustain order, measured with respect to the thermal energy. Naively speaking, it may result quite surprising that angular deviations from perfect alignment as small as those mentioned before can have such big effects. Take $\theta=0.25^{\circ}$ as an example: how can a perturbation of $h_{001}=9|J_{\textit{eff}}|\,\sin(0.25^{\circ})\approx 0.05\,k_{B}T$ create an order parameter as big as $\rho_{S}^{s}\approx 0.25$ for $h/|J_{\textit{eff}}|=9$? We think that the answer to this question relies on the fact that this is not the only energy scale in the system. Indeed, $J_{\textit{eff}}$ is the most relevant energy value here. It guarantees that partial single-charge order is retained within $\beta$ chains at this $T$, so that spins do not respond to the misaligned magnetic field on their own, but as part of finite $\beta$ chains (see, for example, Ref. Sazonov et al., 2012). Something similar happens when dipolar interactions act as a perturbation, explaining that Coulomb forces can sustain single-charge order at big fields when $D/|J_{\textit{eff}}|\approx D/k_{B}T\approx 0.13$.

V Conclusions

Order by thermal disorder is a very peculiar effect that has been studied for almost forty years in the context of frustrated magnetic systems. In spite of this, it has associated very few validated examples in nature. In this work we proposed a scenario (the Ising pyrochlore with antiferromagnetic “all-in–all-out” order in zero field) where it may be detected. An easily-controllable external parameter (a magnetic field in the $[110]$ crystallographic direction) can be used to turn on and off the entropic drive towards order, unmasking thus the effect even in the presence of perturbations coupling to the order parameter. Following a protocol that in principle can be applied to any other uncontrolled energy interactions, we showed that OBD can still be found when dipolar interactions are present, provided that they are weak with respect to the nearest-neighbour exchange interaction ($D<0.13\,|J_{\textit{eff}}|$ for a working temperature $T\approx|J_{\textit{eff}}|/k_{B}$). However, special care needs to be taken in order to ensure a field-misalignment angle from the $[110]$ direction well below $1^{\circ}$. This could be attainable in neutron scattering measurements, one of the techniques we propose to identify order due to thermal fluctuations. Although we know of no ideal material where these experiment could be performed, the measurement of a fast decrease in intensity of the $(220)$ peak in the spin structure factor with increasing field in “all-in–all-out” materials like $\mathrm{Nd_{2}Hf_{2}O_{7}}$ or $\mathrm{Cd_{2}Os_{2}O_{7}}$ would be a clear signature of the OBD mechanism being operative.

Acknowledgements.

We thank H. D. Rosales for help with the structure factor code, M. J. P. Gingras for useful input on the relevant literature, G. Aurelio and M. Kenzelmann for insights about angle control in neutron scattering experiments, and UnCaFiQT (SNCAD) for providing computational resources. The authors acknowledge financial support from Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT) through grants PICT 2013 N*∘2004 and PICT 2014 N∘*2618 and Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) through grant PIP 0446.