Signature of magnetic monopole and Dirac string dynamics in spin ice

Abstract

Magnetic monopoles have eluded experimental detection since their prediction nearly a century ago by Dirac¹. It was recently shown that classical analogues of these enigmatic particles can occur as excitations out of the topological ground state of a model magnetic system, dipolar spin ice². These quasi-particle excitations do not require a modification of Maxwell’s equations, but they do interact through Coulomb’s law and are of magnetic origin. Here, we present an experimentally measurable signature of monopole dynamics. In particular, we show that previous magnetic relaxation measurements in the spin-ice material Dy₂Ti₂O₇ (ref. 3) can be interpreted entirely in terms of the diffusive motion of monopoles in the grand canonical ensemble, constrained by a network of ‘Dirac strings’ filling the quasi-particle vacuum. In a magnetic field, the topology of the network prevents charge flow in the steady state. Nevertheless, we demonstrate the existence of a monopole density gradient near the surface of an open system.

Main

Spin-ice systems^4,5,6, such as Dy₂Ti₂O₇ and Ho₂Ti₂O₇, can be described by a corner-sharing network of tetrahedra forming a pyrochlore lattice of localized magnetic moments, as shown in Fig. 1a. The pairwise interaction is made up of both exchange and dipolar terms

where the rare-earth ions carry a moment of 10 Bohr magnetons, m≈10μ_B, r_{i j} is the distance between spins i and j and S_i is a spin of unit length. The coupling constants are on the 1 K energy scale; for example, for Dy₂Ti₂0₇, |J|m²≈3.72 K and D m²≈1.41 K (ref. 7). These energy scales are 100 times smaller than the crystal field terms⁸ that confine the spins along the axis joining the centres of two adjoining tetrahedra. As a result, on the 1 K energy scale, the moments behave as Ising spins along this axis. Remarkably, within this Ising description, the long-ranged dipolar interactions are almost perfectly screened^7,9 at low temperature, with the result that the low-energy properties are almost identical to those of an effective frustrated nearest-neighbour model with antiferromagnetic interactions¹⁰ of strength J_eff=(5D−J)m²/3. This is equivalent to Pauling’s model for proton disorder in the cubic phase of ice¹¹, which has extensive ground-state entropy and violates the third law of thermodynamics⁶. It successfully reproduces the thermodynamic behaviour of both ice¹² and spin ice¹³ and describes the microscopic properties of the latter to a good approximation. The extensive set of spin-ice states satisfy the Bernal–Fowler ice rules¹⁴: a 3d analogue of the six-vertex model with topological constraint consisting of two spins pointing into and two out of each tetrahedron (2 in–2 out), as shown in Fig. 1a. Flipping one spin breaks the constraint leaving neighbouring tetrahedra with 3 in–1 out and with 3 out–1 in, which constitute a pair of topological defects. Within the nearest-neighbour model, creation of the defect pair costs energy 4 J_eff, whereas further spin flips can move the defects at zero energy cost. It has recently been shown² that including the full dipolar Hamiltonian of equation (1) leads to an effective Coulombic interaction between the topological defects separated by distance r, μ₀q_iq_j/4πr, where μ₀ is the permeability of free space, q_i=±q=±2m/a, and a is the distance between two vertices of the diamond lattice (see Fig. 1); that is, to a Coulomb gas of magnetic monopoles. Standard electromagnetic theory does allow for such excitations¹⁵, which correspond to divergences in the magnetic intensity H, or magnetic moment M, rather than in the magnetic induction: ∇ · B=∇ ·(H+M)=0. On all length scales above the atomic scale, a 3 in–1 out defect seems to be a local sink in the magnetic moment and therefore a source of field lines in H. It can lower its energy by moving in the direction of an external field and therefore carries a positive magnetic charge¹⁶. What is remarkable about spin ice is that it enables the deconfinement of these effective magnetic charges so that they occur in the bulk of the material on all scales, rather than just at the surfaces within a coarse-grained description¹⁵. A two-dimensional equivalent may exist in artificial spin ice, constituting arrays of nanoscale magnets¹⁷.

Figure 1: Spin-ice structure and emergence of monopoles.

a, The magnetic ions (Ho³⁺ or Dy³⁺) lie on the sites of the pyrochlore lattice and are constrained to the bonds of the dual diamond lattice (dashed lines). Local topological excitations 3 in–1 out or 3 out–1 in correspond to magnetic monopoles with positive (blue sphere) or negative (red sphere) charges respectively. b, The diamond lattice provides the skeleton for the network of Dirac strings with the position of the monopole restricted to the vertices. The orientation of the Dirac strings shows the direction of the local field lines in H.

Given the accessibility of these magnetic quasi-particles, the development of an experimental signature is of vital importance and interest. The ‘Stanford’ superconducting coil experiment^2,18 could in principle detect the passage of a single magnetic quasi-particle, but this seems highly unlikely given that the charges have no mass and therefore have diffusive, rather than Newtonian dynamics. A more promising starting point is therefore to look for a monopole signal from magnetic relaxation of a macroscopic sample^3,8,19. The general dynamic behaviour of spin ice is illustrated in Fig. 2 by the magnetic relaxation time, as a function of temperature for Dy₂Ti₂O₇ (ref. 3), taken from bulk susceptibility measurements. The energy scales discussed above give rise to different regimes: the timescale increases in the thermally activated high-temperature regime, entering a quasi-plateau region below 12 K associated with quantum tunnelling processes⁸, before experiencing a sharp upturn below 2 K. The spins are Ising-like below 12 K and the configuration evolves by quantum tunnelling through the crystal field barrier, whereas above this temperature higher crystal field levels are populated and the timescale drops markedly. The quantum tunnelling plateau regime can therefore be well represented by an Ising system with stochastic single spin dynamics and hence should be dominated by the creation and propagation of monopole objects. This is illustrated, in a first approximation, by comparing the data with an Arrhenius law τ=τ₀ exp(2 J_eff/k_BT), as shown by the red curve in Fig. 2. The timescale τ₀ is fixed by fitting to the experimental time at 4 K with J_eff=1.11 K, the value estimated for Dy₂Ti₂0₇ (ref. 7). 2 J_eff is the energy cost of a single, free topological defect in the nearest-neighbour approximation and is half that for a single spin flip. The calculation fits the data over the low-temperature part of the quasi-plateau region, where one expects a significant defect concentration without any double defects (4-in or 4-out), and gives surprisingly good qualitative agreement at lower temperature, as the concentration decreases. Although still in the tunnelling regime, the plateau region corresponds to high temperature for the effective Ising system. Good agreement here provides a stringent test and any theory not fitting must be discarded. The above expression clearly does a good job, enabling us to equate τ₀ with the microscopic tunnelling time. This test therefore already provides very strong evidence for the fractionalization of magnetic charge² and the diffusion of unconfined particles. However, this (or any other) Arrhenius function ultimately fails, underestimating the timescale at very low temperature: although it is possible to fit the data reasonably below 2 K by a single exponential function by varying the barrier height, simultaneous agreement along the plateau and at lower temperature is impossible. The role of the missing Coulomb interaction is therefore clear: although non-confining, it must considerably increase the relaxation timescale by modifying the defect concentration and slowing down diffusion through the creation of locally bound pairs.

Figure 2: Relaxation timescales τ in Dy₂Ti₂O₇: experiment and simulation.

The experimental data (crosses) are from Snyder et al. ³. The Arrhenius law (red line) represents the free diffusion of topological defects in the nearest-neighbour model. The relaxation timescale of the Dirac string network driven by Metropolis dynamics of magnetic monopoles has been obtained for fixed chemical potential (pink filled triangles) and with μ varying slowly to match the defect concentration in dipolar spin ice (blue filled circles). The temperature scale is fixed without any free parameters. Inset: The same data shown in the low-temperature region.

We have tested this idea by directly simulating a Coulomb gas of magnetically charged particles (monopoles), in the grand canonical ensemble, occupying the sites of the diamond lattice. The magnetic charge is taken as q_i=±q. In the grand canonical ensemble, the chemical potential is an independent variable, of which the value in the corresponding magnetic experiment is unknown. In a first series of simulations, we have estimated it numerically by calculating the difference between the Coulomb energy gained by creating a pair of neighbouring magnetic monopoles and that required to produce a pair of topological defects in the dipolar spin-ice model, with parameters taken from ref. 7, giving a configurationally averaged estimate μ/k_B=8.92 K. In a second series of simulations, μ was taken as the value required to reproduce the same defect concentration as in a simulation of dipolar spin ice at temperature T. Here, μ varied only by 3%, with the same mean value as in the first series, showing that our procedure is consistent. The chemical potential used is thus not a free parameter. As the Coulomb interaction is long-ranged, we treat a finite system using the Ewald summation method^20,21. The monopoles hop between nearest-neighbour sites through the Metropolis Monte Carlo algorithm, giving diffusive dynamics, but with a further local constraint: in the spin model a 3 in–1 out topological defect can move at low energy cost by flipping one of the 3-in spins, the direction of the out-spin being barred by an energy barrier of 8 J_eff. An isolated monopole can therefore hop to only 3 out of 4 of its nearest-neighbour sites, dictated by an oriented network of constrained trajectories similar to the ensemble of classical ‘Dirac string’² of overturned dipoles¹⁵. The positively charged monopoles move in one sense along the network, whereas the negative charges move in the opposite direction (see Fig. 1b). The network is dynamically rearranged through the evolution of the monopole configuration. The vacuum for monopoles in spin ice thus has an internal structure: the Dirac strings which, in the absence of monopoles, satisfy the ice rules at each vertex. This structure is manifest in the dynamics and influences the resulting timescales. In fact the characteristic timescale that we compare with experiment comes from the evolution of the network of Dirac strings rather than from the monopoles themselves. Indeed, the monopole autocorrelation time, as extracted from the monopole density–density correlation function (see, for example, ref. 22), turns out to be small for this range of temperature. We locally define the string network by an integer σ=±1, giving the orientation of the Dirac string along each bond of the diamond lattice, and define the autocorrelation function

where t is the Metropolis time and N is the number of bonds (up to N≈25,000). For the initial conditions, we take an ordered network with no monopoles, which we let evolve at temperature T until an equilibrium configuration is attained. This defines t=0. C(t) decays almost exponentially, with characteristic time, τ, that varies with temperature. To avoid initial transient effects, we define τ such that C(τ)=0.8. The time is re-set to zero when C(t) decays beyond 0.01 and the process is repeated many times to give the configurationally averaged decay time. Figure 2 shows a comparison of our simulations with the experimental data of ref. 3. The Metropolis time is again scaled to the experimental time at 4 K and there is again no scale factor on the temperature axis. Data for fixed chemical potential are shown by the pink triangles, whereas data with μ varying are shown by the blue circles. There is a quantitative evolution of the simulation data compared to the nearest-neighbour spin-ice model. Agreement between the experimental and numerical data now looks excellent, showing clearly that the experimental relaxation is due to the creation and proliferation of quasi-particle excitations that resemble classical monopoles in the magnetic intensity H. As the temperature increases, towards the end of the plateau region, a small systematic difference occurs. This is because the spin system can access double defects at finite energy cost, whereas this state corresponds to two like charges superimposed on the same site, which is excluded by the Coulomb interaction. The inset of Fig. 2 shows results at low temperature, illustrating in detail the extent of the improvement in comparison with experiment. Allowing the variation of μ provides a further evolution towards the experimental data, compared with that for fixed μ and the blue circles represent our best numerical results. We now have quantitative agreement between experiment and theory down to low temperature, showing that the Coulomb interactions are responsible for the non-Arrhenius temperature dependence of the relaxation timescales. Differences remain at this level of comparison below 1 K, but to go further would require an even more detailed modelling of spin ice²³ as well as complementary experimental measurements.

Finally we consider the response of monopoles to an external magnetic field, h, placed along one of the [100] directions. Applying such a field to a system for closed-circuit geometry (periodic boundaries), one might expect the development of a monopole current in the steady state¹⁶. This is not the case, at least for the nearest-neighbour model, where we find that a transient current decays rapidly to zero (see Fig. 3a). The passage of a positive charge in the direction of the field reorganizes the network of strings, leaving a wake behind it that can be followed either by a negative charge, or by a positive one moving against the field, with the result that the current stops. This is a dynamic rather than static effect and is not related to confinement of monopole pairs by the background magnetization². Reducing the temperature at finite field, the magnetization saturates around a critical temperature: a vestige of the Kasteleyn transition²⁴, which is unique to topologically constrained systems. Confinement occurs here, as the Zeeman energy outweighs the entropy gain of free monopoles. The transient currents suggest the development of charge separation in an open system. This is indeed the case despite the fact that monopole numbers are not conserved at open boundaries. Figure 3b shows the profile of positive charge density across a sample of size L, with open boundaries, for varying fields. There is a clear build-up of charge over a band of 4–5 lattice spacings, although including long-range interactions may lead to a quantitative change in this value. As the ratio T/h and the monopole density go to zero, the band narrows and the system forms a conventional layer of magnetic surface charge as one expects for any magnetically ordered system¹⁵. In the absence of topological defects, the magnetization is conserved from one layer to another, so that a charge density profile manifests itself as a magnetization profile. The data here suggest charge build-up in a layer several nanometres thick, making it in principle a measurable effect.

Figure 3: Monopole density profile.

a, A magnetic field h is applied at t=0 along the z axis ([100] direction). We show the transient flux of positive charges, Φ, passing through a plane perpendicular to the field, as a function of Metropolis time, t. The simulations are obtained using the nearest-neighbour spin-ice model with periodic boundary conditions (filled squares) and open boundaries, with current measured either at the surface (blue filled triangles) or in the bulk (red filled circles). b, Density of positive defects in the horizontal planes for T=1 K, as a function of z and h (in units of k_BT/m μ₀).