format: text

text:
# Dirty quantum Hall ferromagnets and quantum Hall spin glasses

## I Introduction

Disorder and interactions have competing consequences in quantum Hall ferromagnets (QHFM). In this paper, we study how a disordered impurity potential can give rise to a spin-disordered ground state. We also discuss the influence of disorder on the magnetic and charge response of such a system. We use a classical spin model throughout to describe the quantum Hall ferromagnet.

Coulomb interactions lead to spin correlations in a quantum Hall system. For electrons fully occupying the lowest Landau level (filling fraction $`\nu =1`$), exchange is responsible for a spin-polarised ground state, even in the absence of Zeeman energy. This is the consequence of Hund’s rule as applied to a macroscopically large number of degenerate Landau-level orbitals. The resulting quantum Hall ferromagnet is especially interesting as a system in which the spin configuration and the charge density are closely linked Girvin (1999). At $`\nu =1`$ and if Zeeman energy is large, charge enters the spin-polarised system as minority-spin electrons. However, if Zeeman energy is small or vanishing, the charged excitation of lowest energy is not a bare spin-half electron, but a bound state of an electron with many spin waves. In classical terms, this occurs because the minority spin polarises its local ferromagnetic background, and the composite object may be viewed as a topological excitation, or texture in an ordered ferromagnet — a skyrmionSondhi et al. (1993). Similarly, an antiskyrmion, with topological charge of the opposite sign, is produced when charge is removed from a filled Landau level. In this description, the deviation of local charge density from that of a filled and ferromagnetically polarised Landau level is proportional to the topological density Rajaraman (1989) of the spin configuration.

In a clean system with sufficiently small Zeeman energy, skyrmions or anti-skyrmions can be created at zero temperature on varying the average filling factor from $`\nu =1`$ to larger or smaller values. For a disordered quantum Hall ferromagnet, the coupling of an electrostatic impurity potential to the charge density offers an additional mechanism by which spin textures are nucleated at zero temperature. The physical consequences of this coupling are the subject of this paper, which provides a detailed account of work presented in outline in Ref. Rapsch et al., 2002.

The interplay between disorder and exchange in quantum Hall ferromagnets has been examined previously from several different viewpoints. Fogler and Shklovskii Fogler and Shklovskii (1995) have developed a theory particularly applicable in higher Landau levels. Building on earlier discussions Ando and Uemura (1973), they set out a mean-field treatment in the spirit of Stoner theory. For odd integer filling in the absence of Zeeman coupling, they find a transition between ferromagnetic and paramagnetic ground states with increasing disorder strength. They suggest that this transition should be apparent in transport measurements, in which the ferromagnetic phase is characterised by spin-resolved Shubnikov-de Haas oscillations, and the paramagnet by spin-unresolved oscillations. Experimentally, a transition of this kind is observed with decreasing magnetic field strength Paalanen et al. (1982); Leadley *et al.* (1998), and its sharpness suggests that its origin is indeed cooperative.

Within the Fogler-Shklovskii approach, local moments are all collinear in the ferromagnet and vanish at the transition to the paramagnet. An alternative scenario may arise in the lowest Landau level near $`\nu =1`$, in which the QHFM responds to disorder mainly through the direction rather than the magnitude of its local magnetisation. Some indications that this can happen come from calculations for the fully-polarised ferromagnet at weak disorder. Here, a reduction in spin stiffness with increasing disorder strength has been interpreted by Green Green (1998) as a precursor of a non-collinear phase. Moreover, even weak disorder may nucleate a dilute glass of skyrmions and anti-skyrmions at the maxima and minima of the disordered potential, as discussed by Nederveen and NazarovNederveen and Nazarov (1999) and examined further in the present paper. In addition, Sinova, MacDonald and Girvin Sinova et al. (2000) have shown that, at intermediate disorder strength, both reduced and non-collinear local moments emerge from a numerical solution of Hartree-Fock theory for a model with Coulomb interactions and spatially uncorrelated disorder, while transport properties within a Hartree-Fock theory have been discussed by Murthy Murthy (2001).

In general, we believe that the relative importance of local moment reduction versus the formation of spin textures for dirty QHFM’s will depend on the nature of the disorder. In this work, we concentrate on textures, which are favoured by a smoothly varying impurity potential at $`\nu =1`$. We find that the ferromagnet gives way to a spin glass at strong disorder.

The plan of this paper is as follows. In the next section, we describe a spin model with quenched disorder which is intended to capture the physics of the disordered QHFM. This is followed in section III by a discussion of the phase diagram for the model as a function of filling factor and disorder strength, using heuristic arguments and scaling ideas. The conclusions are supported by the results of a Monte Carlo simulation of a lattice version of the spin model. The technique is outlined in section IV and the results are presented in section V. In sections VI and VII, we discuss the compressibility and conductivity of the system in order to characterise the charge response of spin-disordered ground states.

## II Spin Model

Consider a two-dimensional electron gas in a strong perpendicular magnetic field $`B`$, with Landau-level filling $`\nu `$ close to unity. The electrons are subject to an impurity potential $`V(𝐫)`$ and an electron-electron interaction energy U(r). As a first step, let us omit the exchange interactions and the Zeeman energy. Then the electron density $`\rho (𝐫)`$ is determined by the balance between disorder and interactions, or in other words, screening. We treat this using Thomas-Fermi theory. Such an approximation has been applied by EfrosEfros (1988, 1989) to the comparable problem in spin-polarised Landau levels when $`\nu `$ lies near half-integer values. The ground-state charge density $`\rho (𝐫)`$ at weak disorder is determined by the condition that the Hartree potential should match the chemical potential $`\mu `$ everywhere:

$$\mu =V(𝐫)+U(𝐫𝐫^{})\rho (𝐫^{})d^2𝐫^{}.$$
(1)
This approach is valid in the case where the resulting local filling factor varies smoothly on the scale of the magnetic length, $`l_B=(\mathrm{}/eB)^{1/2}`$, and only has small fractional deviations from $`\nu =1`$, so that

$$\delta \rho (𝐫)\rho (𝐫)(2\pi l_B^2)^1\rho (𝐫).$$
(2)
However, the Thomas-Fermi picture of good local screening does not take in account exchange interactions. Provided electron density fluctuations are small and smoothly varying, ferromagnetic exchange leads locally to a maximal ferromagnetic polarisation of the electron spins. This local magnetisation may vary in space. Denoting its direction by the three-component unit vector $`\stackrel{}{S}(𝐫)`$, spatial fluctuations in spin orientation are linked to electron density by Sondhi et al. (1993); Rajaraman (1989)

$$\delta \rho (𝐫)=\frac{1}{8\pi }ϵ_{ij}\stackrel{}{S}\left(_i\stackrel{}{S}\times _j\stackrel{}{S}\right).$$
(3)
This direct connection is specific to the quantum Hall ferromagnet — a varying electron density implies a variation in the direction of the local magnetisation and *vice versa*. Such spin textures cost exchange energy. So, a proper description of the system must include exchange, impurity and Hartree contributions to the total energy of a dirty quantum Hall ferromagnet. This brings us to study the Hamiltonian

$``$ $`=`$ $`{\displaystyle }({\displaystyle \frac{J}{2}}|\stackrel{}{S}(𝐫)|^2+[V(𝐫)\mu ]\delta \rho (𝐫)`$ (4)
$`+{\displaystyle \frac{U_0}{2}}[\delta \rho (𝐫)]^2)d^2𝐫,`$
where

$$J=\frac{1}{16(2\pi )^{1/2}}\frac{e^2}{4\pi ϵ_0ϵ_rl_B}$$
(5)
is the exchange couplingMoon *et al.* (1995). ($`ϵ_r`$ is the relative permittivity in the semiconductor.) At this point we have chosen for simplicity a short-range Hartree interaction, $`U(𝐫𝐫^{})=U_0\delta (𝐫𝐫^{})`$. We have also absorbed the constant $`U_0/2\pi l_B^2`$ into the chemical potential $`\mu `$.

As mentioned above, we will work with a disordered potential $`V(𝐫)`$ that is smooth on the scale of the magnetic length $`l_B`$. For simplicity, our discussion of this continuum model will assume a Gaussian distribution with standard deviation $`\mathrm{\Delta }`$ and correlation length $`\lambda `$ much larger than the magnetic length $`l_B`$. (Our numerical study will use a lattice model with a bounded distribution with uncorrelated disorder.)

In restricting our study to this model, we neglect quantum fluctuations of $`\stackrel{}{S}(𝐫)`$. This semiclassical approximation is justified for smooth variations, with $`|\stackrel{}{S}(𝐫)|l_\mathrm{B}^1`$. Our aim in the following is to understand the zero-temperature phase diagram of the model defined by Eq. (4), as a function of disorder strength $`\mathrm{\Delta }`$ and average charge density $`\delta \rho `$, the spatial average of $`\delta \rho (𝐫)`$. We will characterise its ground states via their response functions and excitations.

We finish this section by comparing this model with some other examples of disordered systems. As an electron system, it is unusual in that there is an exchange gap for single-particle excitations, even if the ground-state spin configuration $`\stackrel{}{S}(𝐫)`$ does not have long-range ferromagnetic order. This means that the only low-energy excitations involve collective spin modes. As a ferromagnet with quenched disorder, the system is also unusual in several ways. First, the link between spin and charge means that the spin system responds to applied electric fields. We calculate in the following the wavevector-dependent dielectric constant, $`ϵ(q)`$, and compare it with behaviour found in more conventional disordered electron systems. Second, due to the same coupling, spin waves generate an electric dipole moment. This means that spin waves contribute to the optical conductivity $`\sigma (\omega )`$. More generally, the coupling to disorder in this model leaves spin-rotational symmetry intact but breaks time-reversal symmetry. This is in contrast to the effect of random Zeeman fields, which break both symmetries, and to random exchange interactions, which preserve both symmetries.

## III Phase Diagram

We begin with a qualitative discussion of the phase diagram as a function of disorder strength and average charge density. We employ scaling arguments to obtain the phase boundary for breakdown of long-range ferromagnetic order.

The Hamiltonian of Eq. (4) for the continuum model is minimised by the spin configuration which satisfies (appendix A):

$`J^2\stackrel{}{S}\times \stackrel{}{S}`$ $`=`$ $`ϵ_{ij}_iV_\mathrm{H}(𝐫)_j\stackrel{}{S}`$ (6)
$`V_\mathrm{H}(𝐫)`$ $`=`$ $`V(𝐫)\mu +U_0\delta \rho (𝐫)`$ (7)
where $`V_\mathrm{H}`$ the local Hartree potential and $`\delta \rho `$ is defined in terms of the local spin configuration by Eq. (3). This equation for the ground state configuration shows the interplay between exchange and disorder in the model. Since it is difficult to tackle the non-linear equation directly, we proceed using heuristic arguments instead.

The model is characterised by two energy scales: the exchange energy $`J`$ and the disorder strength $`\mathrm{\Delta }`$. There are also two length scales: the correlation length $`\lambda `$ of the disordered potential and

$$L_\mathrm{H}(U_0/J)^{1/2},$$
(8)
which we call the Hartree length. The significance of the Hartree length can be made clear by considering a skyrmion of fixed shape and radius $`R`$ in a clean system. The contribution to its total energy from exchange is $`J`$ and independent of size, while that from Hartree interactions is size dependent, being $`U_0/R^2`$. Comparing these contributions, one sees that exchange dominates on length scales large compared with $`L_\mathrm{H}`$, while Hartree interactions dominate at smaller distances.

Our central hypothesis is that the competition between interaction and disorder in this system is characterised by the Hartree length only. In the following, we also use the limit $`L_\mathrm{H}\lambda `$ as a source of simplifications.

### III.1 Filled Landau level

Let us consider first the effect of disorder on a system at $`\nu =1`$, imposed by setting $`\delta \rho =0`$. Without impurities ($`V(𝐫)=0`$), the system is a perfect ferromagnet. Moreover, there is a threshold Green (1998); Nederveen and Nazarov (1999); Sinova et al. (2000)

$$|V(𝐫)|=4\pi J.$$
(9)
below which an impurity potential is unscreened. It arises because, for any $`\stackrel{}{S}(𝐫)`$, one has $`|\stackrel{}{S}(𝐫)|^28\pi |\delta \rho (𝐫)|`$ and hence

$$\left[4\pi J|\delta \rho (𝐫)|+V(𝐫)\delta \rho (𝐫)\right]d^2𝐫$$
(10)
Thus, if $`|V(𝐫)|`$ is below the threshold $`4\pi J`$ everywhere, the ground state is the perfectly aligned ferromagnet with $`\delta \rho (𝐫)=0`$.

At weak disorder ($`\mathrm{\Delta }J`$), $`|V(𝐫)|`$ in most parts of the system lies below the threshold. With an unbounded potential distribution, the ground state spin configuration therefore consists of a dilute glass of skyrmions and anti-skyrmions, nucleated at rare positions where $`|V(𝐫)|`$ is large, as discussed by Nederveen and NazarovNederveen and Nazarov (1999). Away from these positions we expect that the ferromagnetic order is essentially unaffected by the impurity potential. A careful treatment of this regime is however quite subtle, since the spin deviation due to an isolated skyrmion falls off with distance only as $`r^1`$. We argue in appendix C that long range ferromagnetic order is indeed preserved, and that the internal degrees of freedom of dilute skyrmions and antiskyrmions develop the correlations necessary to ensure this.

In contrast, at strong disorder ($`\mathrm{\Delta }J`$), the charge density provides almost perfect local screening of the disordered potential, so that

$$\delta \rho (𝐫)V(𝐫)/U_0(\mathrm{\Delta }/J1).$$
(11)
Corrections to perfect Thomas-Fermi screening arise at length scales larger than $`L_\mathrm{H}`$, where exchange becomes important. The effect of exchange is to force screening charges to be quantised, since an unquantised charge costs divergent exchange energy in the thermodynamic limit. We can summarise the effect of exchange by dividing the system into regions of area $`L_\mathrm{H}^2`$, finding for each such area the integral

$$Q_{L_H^2}\frac{V(𝐫)}{U_0}d^2𝐫$$
(12)
and adjusting the total screening charge within every region to the integer value closest to $`Q`$. We argue that these integers are predominantly zero in a ferromagnetic phase, and predominantly non-zero in a phase without ferromagnetic order. To see this, consider a well-ordered ferromagnetic phase, in which $`S(𝐫)`$ has small spatial variations around a global direction of magnetisation. In this case, the net topological charge in any region has magnitude much less than one. Conversely, in a phase without such order, unit topological charge will typically accumulate over a region of linear size given by the ferromagnetic correlation length.

This picture leads us to identify the phase boundary of the ferromagnet as the point at which $`Q^2^{1/2}1`$. To estimate $`Q^2`$, note from Eq. (12) that each correlation area of size $`\lambda ^2`$ contributes to $`Q`$ a charge of magnitude $`\lambda ^2\mathrm{\Delta }/U_0`$ and random sign. Over the Hartree area, there are $`L_H^2/\lambda ^2`$ such contributions, so that

$$Q^2^{1/2}(\lambda ^2\mathrm{\Delta }/U_0)(L_H/\lambda )=\lambda L_\mathrm{H}\mathrm{\Delta }/U_0.$$
(13)
The phase boundary is therefore given by

$$\mathrm{\Delta }_cU_0/(\lambda L_\mathrm{H})=J(L_\mathrm{H}/\lambda ).$$
(14)
For $`\mathrm{\Delta }>\mathrm{\Delta }_c`$, the ground state is strongly disordered and has no ferromagnetic order. Within our classical description, the spins are frozen at zero temperature. We can therefore identify this phase as a quantum Hall spin glass.

### III.2 Away from integer filling

Consider next the ground state away from $`\nu =1`$. We examine behaviour at fixed $`\mathrm{\Delta }<\mathrm{\Delta }_c`$, as a function of the average charge density $`\delta \rho \nu 1`$. For $`\delta \rho =0`$ and $`\mathrm{\Delta }<\mathrm{\Delta }_c`$, the system has a net magnetisation, but the introduction of charge, in the form of skyrmions (or anti-skyrmions) disrupts ferromagnetic spin alignment. The size $`R`$ of an isolated skyrmion in a clean system is divergent if Hartree repulsion is not balanced by Zeeman energy. However, the presence of an impurity potential establishes an optimal size, because it is energetically favourable for a skyrmion to locate its charge distribution in randomly occurring potential wells. The typical depth of such a well of radius $`R\lambda `$, averaged over its area, is $`\mathrm{\Delta }(\lambda /R)`$. (The case $`\lambda >R`$ is treated in Ref. Nederveen and Nazarov, 1999.) The Hartree energy is $`U_0/R^2`$. The total energy is hence minimised if

$$RU_0/(\mathrm{\Delta }\lambda )=L_\mathrm{H}(\mathrm{\Delta }_c/\mathrm{\Delta }).$$
(15)
Note that this value of $`R`$ exceeds $`L_\mathrm{H}`$ for $`\mathrm{\Delta }<\mathrm{\Delta }_c`$. Because of this, we expect that exchange energy dominates over the potential energy and that the skyrmion will not be strongly distorted by the random potential. In other words, screening as discussed in the previous subsection does not alter the present argument.

We expect that ferromagnetic order will persist with increasing charge density until such skyrmions overlap. The phase boundary hence lies at

$$\delta \rho _c=L_\mathrm{H}^2(\mathrm{\Delta }/\mathrm{\Delta }_c)^2.$$
(16)
In summary, we have used simple arguments to obtain the phase boundaries between the ferromagnetic and the spin glass as a function of disorder and charge density. The results, Eqs. (14) and (16), are summarised in the schematic phase diagram shown in Fig. 1. We next present results from Monte Carlo simulations which support these predictions.

The discussion we have presented is for a model with a short-range Hartree interaction and for the limit $`l_B\lambda L_\mathrm{H}`$. The central consequence of using a Coulomb form, $`U(r)=e^2/4\pi ϵ_0r`$, in place of a short-range Hartree interaction is to change the length scale derived by comparing Hartree with exchange energies, from $`L_\mathrm{H}(U_0/J)^{1/2}`$ to $`L_\mathrm{H}e^2/4\pi ϵ_0J=16(2\pi )^{1/2}l_B`$. Since in this case $`L_\mathrm{H}`$ is not parametrically larger than $`l_B`$, scaling arguments of the type presented in this section are not justified. In this case, in place of three distinct regimes of disorder strength ($`\mathrm{\Delta }J`$, $`J\mathrm{\Delta }\mathrm{\Delta }_c`$ and $`\mathrm{\Delta }_c\mathrm{\Delta }`$) we expect only two, with $`\mathrm{\Delta }_cJ`$.

## IV Monte Carlo simulations

In order to test the theoretical results derived in the previous section, we use Monte Carlo methods, in combination with simulated annealing, to study a classical Heisenberg model with quenched disorder. In this section, we outline our methodology. Similar techniques have been applied recently to systems without quenched disorder in Ref.Hale et al., 2000. We present our results for the physical response functions of the model in sections V and VI.

We treat a lattice version of Eq. (4) with $`N\times N`$ classical spins $`\stackrel{}{S}_i`$ of unit length on a square lattice, taking periodic boundary conditions. Spins have nearest-neighbour ferromagnetic interactions of strength $`J`$. In addition, the electrical charge $`q_k`$ on each plaquette $`k`$ is calculated from the area on the spin sphere covered by the spins on the corners of that plaquette. This charge density has a local repulsive Hartree interaction of strength $`U_0`$, and is also subject to a uniformly distributed random background potential $`ϵ_k[\mathrm{\Delta },\mathrm{\Delta }]`$. Since we choose this potential independently for each plaquette, the correlation length $`\lambda `$ is set by the lattice spacing.

The Hamiltonian of the lattice model is

$`H`$ $`=`$ $`J{\displaystyle \underset{ij}{}}\stackrel{}{S}_i\stackrel{}{S}_j+{\displaystyle \underset{k}{}}(ϵ_kq_k+{\displaystyle \frac{U_0}{2}}q_k^2)`$ (17)
$`+\gamma (Q_0{\displaystyle \underset{k}{}}q_k)^2,`$
where we have introduced a Lagrange multiplier $`\gamma `$ to bias the system towards a predefined number $`Q_0`$ of charge quanta.

To obtain the ground state spin configuration for a given disorder realisation, we start from a random initial state and anneal using Monte Carlo dynamics and the Metropolis algorithm. After some experimentation, we found the following three-stage protocol to be effective. In the first stage, the temperature is reduced linearly in time from a high temperature $`T_0`$ (several times $`J`$) to $`T_0/10`$, using $`3.5\times 10^5`$ Monte Carlo steps per spin (MCS). In the second stage, the temperature is reduced from $`T_0/10`$ to $`0`$, using $`5\times 10^5`$ MCS. For both these stages, the attempted spin update is an isotropically distributed reorientation. In the third stage, the system is quenched for a further $`5\times 10^4`$ MCS, using as the attempted spin update only small angle reorientations in order to improve the acceptance rate.

We have checked that this algorithm finds the ground state reliably for a weakly disordered system with overall charge neutrality, by doing repeated runs for a given realisation of background potential and using local charge and energy densities to identify states that differ only by a global spin rotation. In strongly disordered systems and those with non-zero average charge density, repeated applications of the algorithm do not reproduce the same state to high precision. Instead, a number of low-lying states are obtained, having a small spread (less than 10%) in their energies and other observables such as their magnetisation. We perform at least five independent simulations for every disorder realisation and pick from the states obtained the one with the lowest energy. For a chosen disorder strength $`\mathrm{\Delta }`$, we generate three or more realisations of the disorder potential and average our results over these realisations. We find that properties such as the ground state energy, the magnetisation and various response functions show only small fluctuations between different disorder realisations for systems of size 32 $`\times `$ 32 spins or larger. Simulations of very much bigger systems are ruled out by constraints on computing time, and we choose $`N=40`$ as the system size for most of our simulations.

The Lagrange multiplier $`\gamma `$ is chosen sufficiently large that a majority of the simulation runs yield a ground state with the desired total charge $`Q_0`$. We do not detect any dependence of our results on the particular value of $`\gamma `$.

## V Magnetic Response

Having outlined our simulation techniques, we next present our results. We start with the magnetic properties of the system. We show that they indicate a transition from the ferromagnet to the spin glass with increasing disorder strength, as expected from our discussions in section III.

### V.1 Magnetisation

We calculate the site-averaged magnetisation

$$M=|\stackrel{}{S}(𝐫)|=N^2|\underset{i}{}\stackrel{}{S}_i|.$$
(18)
Results are shown for a system of 40$`\times `$40 spins with $`U_0=8\pi J`$ in Fig. 2(a). Spins are fully aligned for $`\mathrm{\Delta }/4\pi J<1`$, because our bounded disorder distribution then lies entirely below the threshold, Eq. (9). Increasing $`\mathrm{\Delta }`$ beyond $`4\pi J`$, we obtain a partially polarised ferromagnet. In fact, the magnetisation remains close to its saturated value until $`\mathrm{\Delta }/4\pi J1.8`$, when it starts to drop appreciably. Increasing the disorder strength further, we reach a regime of strong disorder at $`\mathrm{\Delta }/4\pi J3`$, in which there is only a small magnetisation, of the magnitude expected for a spin-glass state of a finite-size system.

These results are consistent with the existence of a phase transition from the collinear, ordered phase at weak $`\mathrm{\Delta }`$ to a spin-glass phase at strong disorder. From the magnetisation curve we pick $`\mathrm{\Delta }_c/4\pi J2.5`$ as our estimate of the critical disorder strength.

### V.2 Susceptibility

We also calculate the uniform susceptibility $`\chi `$ from the response of the ground state to a weak Zeeman field. Our procedure is as follows. First we obtain a ground state without a Zeeman field, using the protocol described in Sec. IV. Then we apply a weak Zeeman field $`\stackrel{}{h}`$ in the direction of the residual magnetisation, by adding to the Hamiltonian the perturbation $`\delta H=\stackrel{}{h}\stackrel{}{M}`$. We find the ground state in the presence of this perturbation by repeating the third stage of the quenching protocol. Taking care to check that our measurements remain within the linear response regime, we extract the susceptibility from the change in the magnetisation $`\delta \stackrel{}{M}`$, using

$$\chi =|\delta \stackrel{}{M}|/|\stackrel{}{h}|.$$
(19)
Results for $`U_0=8\pi J`$ are shown in Fig. 2(b). The susceptibility as a function of disorder strength has a large peak at $`\mathrm{\Delta }/4\pi J=2.5`$. We interpret this as a second indication of a phase transition from the ferromagnetic phase to a disordered phase.

### V.3 Spin Stiffness

We now turn to the spin stiffness $`\rho _s`$ which measures the rigidity of the spin configuration. It is obtained by calculating the energy cost of small amplitude, long-wavelength spin twists in the ground state. More specifically, let the columns of sites in the $`N\times N`$ lattice be labelled by integers $`1,2,\mathrm{}N`$. Then, starting with a ground state obtained as in Sec. IV, we construct a twisted state by rotating all spins on column $`N/2`$ of the system through a small angle $`\theta `$ about an axis $`\stackrel{}{e}`$. Using the third stage of the quenching protocol, we then relax all spins in this twisted state except those on columns $`1`$ and $`N/2`$, which are held fixed. From the difference in energy $`\mathrm{\Delta }E`$ between the initial and final states, we obtain the spin stiffness for rotations about the axis $`\stackrel{}{e}`$, using

$$\rho _s=\mathrm{\Delta }E/2\theta ^2.$$
(20)
Repeating this for different axes of rotation, we calculate the full $`3\times 3`$ symmetric tensor for the spin stiffness, $`\rho _s^{\alpha \beta }`$. As expected, in the ferromagnetically ordered phase one of the principal axes of this tensor lies to a good approximation along the magnetisation direction, and it is convenient in these calculations to choose rotation axes $`\stackrel{}{e}`$ in directions parallel and perpendicular to $`\stackrel{}{M}`$.

The eigenvalues of the spin stiffness tensor are shown as a function of disorder strength in Fig. 2(c). For the fully-polarised ferromagnet, rotations about the magnetisation direction do not alter the spin configuration, and so one eigenvalue is zero while the other two are degenerate, taking the value $`\rho _s=J`$. For the partially-polarised ferromagnet, all three eigenvalues are non-zero, with two remaining degenerate. The spin-glass state, however, has no special spin direction and this magnetic isotropy means that that all three eigenvalues are approximately degenerate. The stiffness is reduced in value from $`J`$ in the fully polarised ferromagnet to approximately $`J/2`$ in the spin glass. (A variational estimate for the stiffness in the spin-glass phase is $`2J/3`$. See appendix A.) The magnetically isotropic phase is observed for $`\mathrm{\Delta }/4\pi J>2.5`$, yielding the same estimate of $`\mathrm{\Delta }_c`$ as our magnetisation and susceptibility data.

### V.4 Spin Correlation Length

The behaviour of spin correlations provides a further way of characterising ground states. In particular, we consider the correlation function

$$C(r)=\frac{1}{N_r}\underset{(r)}{}\stackrel{}{S}_i\stackrel{}{S}_j,$$
(21)
where the sum runs over all $`N_r`$ spin pairs of separation $`r`$. We extract the spin correlation length $`\xi `$ by fitting to the form

$$C^{}(r)=M^2+(1M^2)\mathrm{exp}(r/\xi )$$
(22)
The behaviour of the correlation length as a function of the disorder strength in the spin-glass phase is shown in Fig. 3. From the Harris criterionHarris (1974), we expect it to diverge as $`\xi (\mathrm{\Delta }\mathrm{\Delta }_c)^\nu `$ with $`\nu >2/d=1`$ as the ferromagnetic phase boundary is approached. Our results are consistent with a divergence at $`\mathrm{\Delta }/4\pi J=2.5`$, although it appears that our they are affected by finite-size effects for $`\xi `$ for $`\mathrm{\Delta }/4\pi J<2.7`$. Perhaps because of these finite-size effects, this fit gives a low value for the exponent: $`\nu =0.7`$. Attempts at a similar fit in the ferromagnetic phase are unsuccessful, and indeed the arguments of Appendix C suggest that correlations in this case may decay with a power law.

### V.5 Phase Diagram

Our Monte Carlo results outlined above allow us to distinguish between a ferromagnetic phase and a quantum Hall spin glass. For a system with $`U_0=8\pi J`$, we conclude that a phase transition from a collinear ferromagnet to a spin-disordered phase occurs at $`\mathrm{\Delta }/4\pi J=2.5`$. Repeating these calculations for different disorder strength and charge density, we can map out a phase diagram for the ground state of the system. Results for $`U_0=4\pi J`$ are shown in Fig. 4.

Comparing this with the phase diagram (Fig. 1) predicted from our heuristic arguments, we see that the two are qualitatively very similar, even though the simulations are carried out for $`(L_H/\lambda )=\sqrt{4\pi }`$ while our scaling picture applies in the limit $`(L_H/\lambda )1`$. An idea of the dependence on $`L_H/\lambda `$ is given by contrasting results at the two values of $`U_0`$ studied. The critical disorder $`\mathrm{\Delta }_c`$ at which the neutral system loses ferromagnetic order is $`\mathrm{\Delta }_c/4\pi J2.5`$ for $`U_0=8\pi J`$, and reduces to $`\mathrm{\Delta }_c/4\pi J2.2`$ for $`U_0=4\pi J`$, in qualitative agreement with the scaling behaviour we expect from Eq. (14).

It is interesting to note that there is a density range over which disorder may *stabilise* the ferromagnet, by limiting the size of the nucleated spin textures. This range is however very narrow: $`|\delta \rho |10^2`$ in units of charge per plaquette.

We mention in passing that we have not searched extensively for a skyrmion crystal, expected at a finite charge density in the weak disorder limit but presumably unstable for $`\mathrm{\Delta }0`$.

## VI Dielectric response

We have so far discussed the phase diagram of the system in terms of its magnetic correlations and response. We next study its charge response.

We examine first the dielectric response of the partially polarised ferromagnet and spin glass, characterised at zero frequency by the wavevector-dependent dielectric constant $`ϵ(q)`$ or by the compressibility $`\kappa (q)`$, related via

$$\kappa (q)=q^2ϵ(q)ϵ_0/e^2,$$
(23)
where $`e`$ is the electron charge and $`ϵ_0`$ is the permittivity of the medium. More precisely, we apply a periodic modulation, $`V(𝐫)V(𝐫)+V_1\mathrm{cos}(𝐪𝐫)`$, to the potential in Eq. (4). As a result, the ground-state electron density changes according to $`\delta \rho (𝐫)\delta \rho (𝐫)+\delta \rho _1(𝐫)`$. Since the system is disordered, $`\delta \rho _1(𝐫)`$ contains many Fourier components, but after averaging the linear response is

$$\delta \rho _1(𝐫)=V_1\kappa (q)\mathrm{cos}(𝐪𝐫),$$
(24)
which constitutes our definition of $`\kappa (q)`$.

Our numerical results for $`\kappa (q)`$ are displayed in Fig. 5, where we study systems deep in the spin-glass phase, setting $`\mathrm{\Delta }/4\pi J=3`$ and $`\delta \rho (𝐫)=0`$. We compare behaviours at $`U_0/4\pi J=1`$ and at $`U_0/4\pi J=2`$, in each case combining data from lattices of size $`40^2`$ and $`56^2`$ in order to maximise wavevector resolution.

We find that, at small $`q`$, $`\kappa (q)`$ is independent of $`U_0`$ and quadratic in $`q`$. At large $`q`$, it is independent of $`q`$, varying roughly as $`U_0^1`$. We can understand these results for $`\kappa (q)`$ using the approach we employed to discuss the phase diagram. The Hartree length $`L_\mathrm{H}`$ again plays an important role.

For $`qL_\mathrm{H}^1`$, exchange may be neglected and we see from Eq. (12) that

$$\kappa (q)U_0^1(qL_H1),$$
(25)
in agreement with the $`q`$-independent value obtained for $`\kappa (q)`$ at large wavevectors from our simulations, and with the $`U_0`$ dependence of these values.

Exchange becomes important at longer wavelengths $`qL_\mathrm{H}^1`$. To estimate $`\delta \rho _1(𝐫)`$ for small $`q`$, we suppose that it arises primarily from spin rotations which have amplitude $`\delta \theta `$ and wavevector $`𝐪`$. To be specific, consider for the periodic perturbation $`V_1\mathrm{cos}(qx)`$ a region of size $`2L\times L`$ with $`L=\pi /q`$. The perturbation causes a net movement of charge $`\delta Q_1`$ in the $`x`$ direction, from one side of this region to the other. We choose to focus on those spins lying on the line of length $`L`$ which is parallel to the $`y`$ axis and which divides the region into two equal halves of size $`L\times L`$. (See Fig. 6.) We suppose, in a variational spirit, that the charge movement induced by the perturbation involves simply a rigid rotation of these spins, through an angle $`\delta \theta `$. With this assumption, let us estimate the charge transfer generated by such a rotation, and the exchange energy that it costs. To do so, it is useful to regard the spin configuration as a map of the line we have defined onto a trajectory on the spin sphere. For a disordered spin configuration, this trajectory resembles a random walk. The end points of the walk are randomly placed on the spin sphere if $`L`$ is large compared to the spin correlation length $`\xi `$, and the distance $`D`$ between them on the surface of the spin sphere is therefore typically $`𝒪(1)`$. Under a rigid rotation of spins on the line, the trajectory is displaced rigidly around the spin sphere. The charge transferred across the real-space line is proportional to the area swept out on the spin sphere by the trajectory during this displacement. We therefore arrive at the estimate

$$\delta Q_1D\delta \theta \delta \theta .$$
(26)
Choosing the phase and axis of rotation appropriately in each such region, the change in potential energy density arising from the rotation is $`V_1\delta Q_1/L^2=V_1q^2\delta \theta `$ while the associated change in exchange energy density is $`Jq^2\delta \theta ^2`$. Choosing $`\delta \theta `$ to minimise the total energy, we find $`|\delta \rho _1(𝐫)|V_1q^2/J`$ and hence

$$\kappa (q)q^2/J=U_0^1(qL_\mathrm{H})^2(qL_H1)$$
(27)
This conclusion is again consistent with our numerical results.

Alternatively, we can arrive at this form for the compressibility from scaling considerations. In general, we may expect that the compressibility is described by the scaling form

$$\kappa (q)=U_0^1f(qL_H),$$
(28)
where $`f(x)`$ approaches a constant at large $`x`$. Our central hypothesis is that exchange dominates at small wavevectors, and so $`\kappa (q)`$ should be independent of $`U_0`$ as $`q0`$. This implies that $`f(x)x^2`$ for $`x1`$, leading to the form given in Eq. (27) for the compressibility at small $`q`$.

To summarise, the system has a metallic response to a perturbing potential at large wavevectors, with $`\kappa (q)`$ independent of $`q`$. However, at small wavevectors, it behaves like an insulator, with $`ϵ(q)`$ independent of q.

We note that, for an infinite system at exactly $`q=0`$ and zero frequency, $`\kappa (0)`$ should be proportional to the thermodynamic density of states. So, except in the fully-polarised ferromagnet, we expect that $`\kappa (0)`$ remains finite as we take the temperature to zero. To reconcile this general expectation with our results, one must remember that for finite-size systems considered here at zero temperature, the discreteness of charge implies that almost all disorder realisations have no response to an infinitesimal change in the chemical potential, while there is a divergent response from those realisations for which ground states from two different charge sectors are degenerate.

## VII Optical Conductivity

Finally, we consider the optical conductivity $`\sigma (\omega )`$ at frequency $`\omega `$. Within our treatment, spin waves are the only excitations that contribute. Since they are not topological excitations and so do not carry net electrical charge, the dissipative conductivity vanishes in the low-frequency limit. Spin waves do, however, give rise to local charge fluctuations and a fluctuating electric dipole moment which couples to an oscillating external electric field, generating dissipation at finite frequency.

For a fully-polarised quantum Hall ferromagnet, GreenGreen (1999) has shown that the spin-wave contribution to optical conductivity is very small. In the non-collinear quantum Hall spin glass, the contribution may be larger due to the presence of a finite charge density in the disordered ground state. Furthermore, the low-energy dynamics of the collinear quantum Hall ferromagnet is qualitatively different from that of a quantum Hall spin glass, since while spin waves in a collinear background have a quadratic energy dispersion, those in a non-collinear background have a linear dispersionHalperin and Saslow (1977), $`\omega =cq`$ at small $`q`$, with velocity

$$c=(\rho _s/\chi )^{1/2}.$$
(29)
Of the three polarisation modes, one is expected to remain gapless even in the presence of a Zeeman coupling.

We now calculate the spin-wave contribution to the optical conductivity for a ground state with non-collinear spins, and then estimate its magnitude in realistic systems.

The rotation of spins from their ground-state orientation in the presence of a spin wave may be parameterised by a vector $`\stackrel{}{p}(𝐫)`$, with at first order

$$\stackrel{}{S}\stackrel{}{S}\stackrel{}{p}(𝐫)\times \stackrel{}{S}.$$
(30)
Using Eq. (3), we show in appendix A that this induces a change in the electron density

$$\delta \rho (𝐫)\delta \rho (𝐫)\frac{ϵ_{ij}}{4\pi }_i\stackrel{}{p}_j\stackrel{}{S}.$$
(31)
This linear dependence on $`\stackrel{}{p}`$ should be contrasted with spin waves in a collinear ferromagnet, for which the density change is $`O(|\stackrel{}{p}|^2)`$.

To estimate the contribution to the conductivity, we note that a spin wave will couple to an external electric field through the net dipole moment $`𝐏`$ that it induces. The coupling appears in the Hamiltonian in the form

$$\delta H=P_x_x$$
(32)
for an electric field $`_x`$ in the $`x`$-direction. From Fermi’s golden rule, the power $`\mathrm{\Pi }(\omega )`$ absorbed from an oscillating electric field of frequency $`\omega `$ is

$`\mathrm{\Pi }(\omega )`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{}^2}}_x^2|f|P_x|0|^2\mathrm{}\omega g(\omega )`$ (33)
where $`|0`$ is the ground state and $`|f`$ is a state with a single spin wave excited, of excitation energy $`\mathrm{}\omega `$, and $`g(\omega )`$ is the spin wave density of states in frequency. From Eq. (29), we have

$$g(\omega )=\frac{3\omega }{2\pi c^2}L_xL_y,$$
(34)
including three polarisations for a system with linear dimensions $`L_x`$ and $`L_y`$.

We estimate for the matrix element appearing in Eq. (33) by considering the dipole moment induced by a spin wave. We start our discussion using the form taken by a spin wave mode in the absence of ground-state disorder, $`\stackrel{}{p}(𝐫)=\alpha \stackrel{}{p}_0\mathrm{cos}(𝐪𝐫)`$, where $`\stackrel{}{p}_0`$ is a unit vector defining the axis of spin rotations and $`\alpha `$ specifies the amplitude. The charge density induced by an excitation of this type in a disordered ground state with spin correlation length $`\xi `$ has, from Eq. (31), magnitude $`\alpha q/\xi `$. We expect this to fluctuate with a random sign over the length scale $`\xi `$. The spin wave therefore induces electric dipoles of magnitude $`e\xi (q\xi \alpha )`$ in each correlation area $`\xi ^2`$. Averaging over a $`L_x\times L_y`$ system, we find

$$P_x^2\left(e\xi ^2q\alpha \right)^2\frac{L_xL_y}{\xi ^2}e^2\alpha ^2\left(\frac{\omega \xi }{c}\right)^2L_xL_y.$$
(35)
where we have substituted $`q=\omega /c`$.

It remains to determine the amplitude $`\alpha `$ for a single quantum excitation. We show in appendix B that

$$\alpha ^2=\mathrm{}/\chi \omega L_xL_y.$$
(36)
Combining factors and dropping numerical coefficients, the absorbed power is

$$\mathrm{\Pi }(\omega )_x^2\frac{e^2\xi ^2\chi }{\rho _s^2}\omega ^3L_xL_y.$$
(37)
This is an ohmic response, $`\mathrm{\Pi }/L_xL_y=\sigma _x^2/2`$, with conductivity

$$\sigma (\omega )\frac{e^2}{h}\left(\frac{\omega \xi }{c}\right)^2\frac{\mathrm{}\omega }{\rho _s}$$
(38)
Let us estimate the magnitude of this spin wave conductivity. In the spin glass phase, we use our numerical results to estimate $`\rho _sJ/2`$. Also, the numerical results for the lattice shows that $`\chi _{\mathrm{lattice}}`$ of the order of unity in units of $`l_B^2/J`$. The continuum magnetisation is related to the lattice spins by: $`\stackrel{}{m}\stackrel{}{S}\mathrm{}/l_B^2`$. This means that the continuum susceptibility is $`\chi \mathrm{}^2/Jl_B^2`$. Combining these factors, we have

$$\sigma (\omega )\frac{e^2}{h}\left(\frac{\xi }{l_B}\right)^2\left(\frac{\mathrm{}\omega }{\rho _s}\right)^3.$$
(39)
Taking $`\xi =10l_B`$ and $`J=`$ 4-8 kelvin (correcting our earlier valueRapsch et al. (2002)), we find for a frequency of 1 GHz, $`\sigma (\omega )`$ ($`10^3`$-$`10^4`$) $`e^2/h`$. Unfortunately, variable-range hoppingPolyakov and Shklovskii (1993) seems likely to mask this contribution to $`\sigma (\omega )`$.

## VIII Summary and discussion

We have investigated the competition between exchange interactions and disorder in quantum Hall ferromagnets at or near integer filling and at zero temperature. Our approach is tailored to the limit of a smoothly-varying impurity potential.

We find that the ferromagnetic state is destroyed by strong disorder through the creation of skyrmion/anti-skyrmion pairs, or by a finite density of either skyrmions or anti-skyrmions at filling factors sufficiently far from $`\nu =1`$. This behaviour, anticipated from simple scaling arguments, is confirmed in numerical studies of ground state spin configurations, obtained from slow Monte Carlo cooling of an initially random high temperature phase. The disordered phase is identified as a quantum Hall spin glass by the absence of long-range order, the presence of non-vanishing local magnetic moments and a finite spin stiffness.

The quantum Hall spin glass has a zero-frequency dielectric response which interpolates between that of an insulator at small wavevectors, and that of a metal at large wavevectors. It supports gapless spin wave modes that couple to electric fields through a finite dipole moment and contribute to the optical conductivity of the system.

Possible experimental signatures of the phenomena we have discussed follow both from the behaviour of the magnetisation and from the nature of excitations. Measurements of the Knight shift in nuclear magnetic resonance (NMR)Barrett et al. (1995); Khandelwal et al. (2001) provide information on the distribution, sampled in space, of the spin polarisation component parallel to the applied magnetic field, while polarisation-resolved absorption spectroscopy Aifer et al. (1996); Zhitomirsky et al. (2002) can be used to determine the average spin polarisation. Impurity effects are most characteristic at $`\nu =1`$, where they result in a reduced spin polarisation, reaching zero in the spin glass at zero Zeeman coupling, and a finite width to the polarisation distribution. We note that spin polarisation which remains unsaturated even at $`\nu =1`$ is found in absorption spectroscopy,Aifer et al. (1996); Zhitomirsky et al. (2002) and that a broad Knight shift distribution is measured in low-temperature NMR.)Khandelwal et al. (2001) In addition, the existence of excitations at energies lower than the Zeeman gap is characteristic of the partially polarised ferromagnet, and has been advanced as an explanation for a coupling observed in NMR experiments between radio-frequency magnetic fields and the electron system.Kuzma et al. (1998) On the other hand, at filling factors away from $`\nu =1`$, neither local probes nor the nature of excitations distinguish sharply between a skyrme crystal in a clean system and the partially polarised ferromagnet or spin glass induced by disorder.

We are grateful for discussions with N. R. Cooper and S. L. Sondhi. The work was supported in part by the EPSRC under Grant GR/J78327 (JTC), and by the Royal Society (DKKL).

## Appendix A Spin waves and charge fluctuations

Consider a small rotation $`𝖱(\stackrel{}{p})`$ of the spins around $`\stackrel{}{p}`$, so that spin directions transform according to $`\stackrel{}{S}𝖱\stackrel{}{S}`$, with

$$𝖱=e^𝖯\mathrm{and}[𝖯]^{ab}=ϵ^{abc}p^c.$$
(40)
To second order, $`\delta \stackrel{}{S}\stackrel{}{p}\times \stackrel{}{S}+\stackrel{}{p}\times (\stackrel{}{p}\times \stackrel{}{S})/2`$. The charge density deviation from $`\nu =1`$ is given by $`\delta \rho _0=ϵ_{ij}\stackrel{}{S}(_i\stackrel{}{S}\times _j\stackrel{}{S})/8\pi `$. In the presence of an additional spin rotation, it becomes

$$\delta \rho _0+\delta \rho _1=\frac{ϵ_{ij}}{8\pi }𝖱\stackrel{}{S}\left[_i(𝖱\stackrel{}{S})\times _j(𝖱\stackrel{}{S})\right]$$
(41)
so that

$$\delta \rho _1=\frac{ϵ_{ij}}{4\pi }_i\stackrel{}{p}_j\stackrel{}{S}+\frac{ϵ_{ij}}{8\pi }_i\stackrel{}{p}_j(\stackrel{}{p}\times \stackrel{}{S})+O(p^3),$$
(42)
obtained using $`(𝖱^1_i𝖱)\stackrel{}{S}(_i𝖯[𝖯,_i𝖯]/2)\stackrel{}{S}`$ $`=`$ $`\stackrel{}{S}\times _i\stackrel{}{p}+`$ $`(_i\stackrel{}{p}\times \stackrel{}{p})\times \stackrel{}{S}/2`$. We can drop higher-order terms if $`|\stackrel{}{p}|1`$ and $`|\stackrel{}{p}||\stackrel{}{S}|`$, so that $`\delta \rho _1/\delta \rho _01`$.

To first order in $`p`$, the component of $`\stackrel{}{p}`$ parallel to $`\stackrel{}{S}`$ does not affect the charge density. At this order, the charge fluctuation $`\delta \rho _1`$ can be written as

$$\delta \rho _1=ϵ_{ij}_i\stackrel{}{p}_j\stackrel{}{S}.$$
(43)
From continuity, $`\delta \dot{\rho }_1+_iJ_i=0`$, the current density is

$$J_i=ϵ_{ij}\dot{\stackrel{}{p}}_j\stackrel{}{S}\text{+ divergence-free part}$$
(44)
Note that we have only identified the transport current, and not any circulating current in the bulk.

Let us now consider the energy cost of spin rotations for the Hamiltonian in Eq. (4). The change in the exchange energy density $`H_J`$ is

$`\delta H_J`$ $`=`$ $`J_i\stackrel{}{p}(_i\stackrel{}{S}\times \stackrel{}{S})+\delta H_J^{(2)},`$
$`\delta H_J^{(2)}`$ $`=`$ $`{\displaystyle \frac{J}{2}}[|\stackrel{}{p}|^22(\stackrel{}{p}\stackrel{}{S})^2+(_i\stackrel{}{p}\stackrel{}{S})_i(\stackrel{}{p}\stackrel{}{S})`$ (45)
$`(\stackrel{}{p}\stackrel{}{S})(_i\stackrel{}{p}_i\stackrel{}{S})].`$
The change in potential energy density is

$$\delta H_\rho =V_\mathrm{H}(𝐫)\delta \rho _1(𝐫)+\frac{1}{2}\delta \rho _1(𝐫)U(𝐫𝐫^{})\delta \rho _1(𝐫^{})d^2𝐫^{}$$
(46)
where

$$V_\mathrm{H}(𝐫)=V(𝐫)\mu +U(𝐫𝐫^{})\delta \rho _0(𝐫^{})d^2𝐫^{}$$
(47)
is the local Hartree potential and $`\mu `$ is the chemical potential.

An equation satisfied by the spin configuration in the ground state is obtained from requiring that the total energy is unaffected by the rotation $`\stackrel{}{p}`$ to first order. This means that we have to balance the first-order terms in the expressions for $`\delta H_J`$ and $`\delta H_\rho `$, giving

$$_i\left[J_i\stackrel{}{S}\times \stackrel{}{S}ϵ_{ij}V_\mathrm{H}(𝐫)_j\stackrel{}{S}\right]=0$$
(48)
for the spins in the ground state. This non-linear equation demonstrates the competition between ferromagnetic exchange and Thomas-Fermi screening in a disordered QHFM.

From the second-order contributions $`\delta H_J^{(2)}`$ in the exchange energy $`\delta H_J`$, we can make a variational estimate of the spin stiffness in the spin-disordered phase. Suppose that $`\stackrel{}{p}`$ and $`\stackrel{}{S}`$ are uncorrelated. Then the second-order terms average to

$`(\stackrel{}{p}\stackrel{}{S})^2=(_i\stackrel{}{p}\stackrel{}{S})_i(\stackrel{}{p}\stackrel{}{S})`$ $`=`$ $`{\displaystyle \frac{1}{3}}|\stackrel{}{p}|^2`$
$`(\stackrel{}{p}\stackrel{}{S})(_i\stackrel{}{p}_i\stackrel{}{S})`$ $`=`$ $`0.`$ (49)
We can define a disorder-averaged spin stiffness from $`\delta H_J^{(2)}=\stackrel{~}{\rho }_s|\stackrel{}{p}|^2/2`$ so that

$$\stackrel{~}{\rho }_s=2J/3.$$
(50)
## Appendix B Long wavelength spin waves

In this appendix, we review the results of Halperin and SaslowHalperin and Saslow (1977) for hydrodynamic excitations, and of GinzburgGinzburg (1978) for elementary excitations, and adapt these for our purposes.

Both theories deal with long-wavelength, low-energy excitations of a disordered spin system. The microscopic magnetisation density and energy density can be coarse-grained over areas $`A`$ to mean values; in the case of the magnetisation one has

$$\stackrel{}{m}(𝐫)=\frac{1}{A}\underset{iA}{}\stackrel{}{S}_i.$$
(51)
Within a hydrodynamic theory, the low-energy dynamics of the system is determined entirely these coarse-grained quantities, which obey conservation laws and are assumed to fluctuate slowly on the timescale set by local relaxation rates.

Hydrodynamic spin fluctuations have equations of motion which may be written in terms of the rotation angle $`\stackrel{}{p}`$, as used in in Eq. (30), and the local magnetisation, with the form

$`{\displaystyle \frac{p_\alpha }{t}}=\chi ^1m_\alpha `$ (52)
$`{\displaystyle \frac{m_\alpha }{t}}=\rho _s^2p_\alpha `$ (53)
where $`\rho _s`$ is the spin stiffness and $`\chi `$ is the uniform magnetic susceptibility. Using the same variables, the free energy of the system is

$`\mathrm{\Delta }F(\stackrel{}{m},\stackrel{}{p})`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha =1}{\overset{3}{}}}{\displaystyle d^2r\left(\frac{1}{\chi }m_\alpha ^2+\rho _s|p_\alpha |^2\right)}`$ (54)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha =1}{\overset{3}{}}}{\displaystyle d^2r\left(\chi \dot{p}_\alpha ^2+\rho _s|p_\alpha |^2\right)}.`$
A long-wavelength treatment of elementary excitations leads to the same equations of motion and an equivalent expression for the energy of the system, in which $`\rho _s`$ and $`\chi `$ are ground-state quantities. We have argued in section V that both are finite at zero temperature in the spin-glass phase. These equations are reminiscent of the dynamics of the spin waves in a Heisenberg antiferromagnet. They lead to spin waves with the dispersion relation

$$\omega =cq;c=(\rho _s/\chi )^{1/2}.$$
(55)
Considering a single spin-wave mode, with

$$\stackrel{}{p}=\alpha \stackrel{}{p}_0\mathrm{cos}(𝐪𝐫\omega t)(|\stackrel{}{p}_0|=1),$$
(56)
we obtain the amplitude $`\alpha `$ for one quantum by setting the energy to $`\mathrm{}\omega `$, arriving at Eq. (36).

## Appendix C Dilute Skyrmions

In this appendix, we study in more detail the interactions between a dilute set of skyrmions at weak disorder. We show how the relative orientations of internal degrees of freedom of the skyrmions are determined by the energetics of the system.

For a clean system without Hartree interactions, the Hamiltonian reduces to the O(3) non-linear sigma model. The skyrmions are of the Belavin-Polykov type. They do not interact. Hartree interactions alone will lead to a divergent skyrmion size in the absence of a Zeeman field. This is prevented in the presence of a disordered potential, as pointed out in section III. However, the presence of an inhomogeneous potential also means that the skyrmions interact. We discuss here whether a physical picture of isolated skyrmions at weak disorder is justified. In particular, since the spin deviation due to a skyrmion falls off as $`1/r`$ with distance $`r`$, one should ask whether the ferromagnetic polarisation is strongly affected by a collection of dilute skyrmions.

For the purposes of this appendix, we choose a simple disorder distribution in which pinned skyrmions and anti-skyrmions are nucleated by isolated wells and barriers of circular shape, positioned randomly on the plane. The depth (or height) of these potentials exceeds the threshold, Eq. (9), so that they nucleate skyrmions (or anti-skyrmions). This model disorder distribution has the advantage that we can use the Belavin-Polyakov solutions of the non-linear sigma model as a starting point for our analysis.

Consider first the case in which the chemical potential is sufficiently high that there are no anti-skyrmions. A set of Belavin-Polyakov skyrmions is described by an analytic function $`w(z)`$ of the position $`z=x+iy`$ in the plane

$$w(z)=\underset{i=1}{\overset{n}{}}\frac{za_i}{zb_i},$$
(57)
where $`n`$ is the number of skyrmions. The function $`w(z)`$ parametrises the spin configuration via

$`S_x+iS_y`$ $`=`$ $`2w/(1+|w|^2)`$
$`S_z`$ $`=`$ $`(1|w|^2)/(1+|w|^2)`$ (58)
with the boundary condition $`S_x=1`$ at infinity. For a single skyrmion ($`n=1`$), the density profile is

$$\delta \rho _{n=1}(z)\frac{|_zw|^2}{\pi (1+|w|^2)^2}=\frac{d_i^2/4\pi }{(|zz_i|^2+d_i^2/4)^2}$$
(59)
where $`z_i=(a_i+b_i)/2`$ can be interpreted as the position of the skyrmion and $`d_i=|a_ib_i|`$ as its diameter. There remains an internal phase degree of freedom $`\varphi _i`$, defined from $`(a_ib_i)=d_i\mathrm{exp}(i\varphi _i)`$. In multiple skyrmions configurations the relative phases are important.

In a dilute glass, skyrmions are typically separated by distances large compared to their size, and the position and diameter of each skyrmion is determined separately by the balance of Hartree and potential energy in the vicinity of the potential well on which it is centered. In the region far from any well, the Belavin-Polyakov form is a good description for the spin configuration, because it minimises exchange energy. The parameters $`z_i`$ and $`d_i`$ are fixed by the spin configuration within the $`i`$th well, while the values of $`\varphi _i`$ remain to be determined. The density profile in this region is

$$\delta \rho ^{}(z)=\frac{|w|^2}{\pi (1+|w|^2)^2}\left|\underset{i}{}\frac{a_ib_i}{(za_i)(zb_i)}\right|^2$$
(60)
where the prime indicates the restricted area in which this form applies. The interference terms in the squared sum make a contribution to the total charge in the region far from the wells which amounts to a pairwise coupling between the internal degrees of freedom of different skyrmions. The integrated density is

$`I`$ $`=`$ $`{\displaystyle ^{}}{\displaystyle \frac{d^2z}{\pi }}{\displaystyle \frac{|w|^2}{(1+|w|^2)^2}}\times `$ (61)
$`\text{Re}{\displaystyle \underset{ij}{}}{\displaystyle \frac{(a_ib_i)^{}(a_jb_j)}{(z^{}a_i^{})(z^{}b_i^{})(za_j)(zb_j)}}`$
In the dilute limit, we can take $`a_ib_i=z_i`$. We can also approximate $`|w|^21`$ away from the skyrmions. Then

$`I`$ $``$ $`\text{Re}{\displaystyle ^{}}{\displaystyle \frac{d^2z}{4\pi }}{\displaystyle \underset{ij}{}}{\displaystyle \frac{d_id_je^{i(\varphi _i\varphi _j)}}{(z^{}z_i^{})^2(zz_j)^2}}`$ (62)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i>j}{}}{\displaystyle \frac{d_id_j\mathrm{cos}(\varphi _i\varphi _j)}{|z_iz_j|^2}}`$
This contribution to the charge in the region between potential wells adds to the energy of the system, since the electrostatic potential here does not exceed threshold \[Eq. (9)\]. It represents charge density that is removed from the cores of skyrmions, as a result of overlap with the tails of other distant skyrmions. We see that in order to minimise the energy we should choose the skyrmion phases $`\varphi _i`$ so as to minimise $`I`$. This is done by arranging $`\varphi _i`$ so that the variable $`d_ie^{i\varphi _i}`$ sums locally to small values. With such correlations, we expect that long-range ferromagnetic order survives in the presence of dilute skyrmions.

A similar treatment of a mixed system containing both skyrmions and antiskyrmions leads to interactions of the same kind amongst the skyrmions, and separately amongst the antiskyrmions, but without coupling between the two species at leading order in inverse density.


text:
# Magnetization reversal times in the 2D Ising model

## I Introduction

Activated processes that can be described with some type of Brownian dynamics are abundant in the world around us. Well-known examples are the nucleation of droplets in an undercooled gas or of crystals in an undercooled liquid, chemical reactions and the escape of a protein from a misfolded state. A prototype system to study such phenomena numerically is the well-known Ising model. Above the so-called critical temperature, in absence of an external magnetic field, up- and down-pointing spins are roughly equally abundant. Below the critical temperature, the system prefers to be in either of two states: one state with a positive magnetization in which most spins are pointing up, and one state with a negative magnetization. As long as the system size remains finite, reversals of the magnetization - transitions between positive and negative magnetization - are possible and will occur at a certain average frequency. These processes are activated, since configurations with magnetization close to zero have a higher free energy than typical configurations with a magnetization close to either of the equilibrium values.

In this manuscript, we study the time scales associated with magnetization reversal. A theoretical framework is outlined which is generally applicable to activated processes in systems with Brownian type dynamics, and compared to high-accuracy computer simulations. From a practical point of view magnetization reversals are also of great interest because of applications in memory devices and the like. One wants to have rapid switching of magnetization under reversals of an external field, but no spontaneus reversals of the magnetization, even if the external field has been turned off. In the literature much attention has been paid to reversal time distributions in the presence of a driving field , but spontaneous reversals in the absence of a field have hardly been studied. Here we consider the latter case for the prototypical case of an Ising model with periodic boundary conditions. We identify the leading scenario for reversals of the magnetization and show that the process may be described to a good approximation by a one-dimensional diffusion process over a potential barrier.

Our manuscript is organized as follows. In section II, we describe the model that we study in detail. Next in section III, we outline the theoretical framework which is generally applicable to activated processes in systems with Brownian type dynamics. We then apply this framework to our prototypical model —magnetization reversal in the Ising model. In section IV we compare the theoretical predictions with high-accuracy computational results.

## II Detailed description of the model

We consider the Ising model on a $`B\times L`$ rectangular lattice with periodic (helical) boundary conditions, with the Hamiltonian

$$H=J\underset{i,j}{}\sigma _i\sigma _j,$$
(1)
in which $`\sigma _i=\pm 1`$ is the spin at site $`i`$ and $`J`$ is the coupling constant. The summation runs over all pairs of nearest-neighbor sites; those of site $`i`$ are $`j=i\pm 1`$ modulo $`N`$ and $`j=i\pm B`$ modulo $`N`$, with $`N=BL`$. The magnetization is defined as $`M_i\sigma _i`$; it can take values $`M=N,N+2,\mathrm{},N`$; all through this manuscript, we restrict ourselves to systems in which both $`B`$ and $`L`$ are even. As a consequence, $`M`$ takes only even values, and summations over a range of possible magnetizations only run over even numbers, with an increment of 2.

A well-known property of the model is that at sufficiently low temperatures, the distribution of the magnetization $`M_i\sigma _i`$ becomes bimodal: the spins in a configuration tend to align around the two preferred values $`\pm M_0`$. For infinite systems this occurs at all temperatures below the critical temperature $`T_c=J/(0.44069k_B)`$, for finite systems this range starts at slightly different temperatures.

The system evolves in time according to single-spin-flip dynamics with Metropolis acceptance probabilities . If $`C_i`$ is the configuration after $`i`$ proposed spin flips, a trial configuration $`C_{i+1}^{}`$ is generated by flipping a single spin at a random site. This trial configuration is then either accepted ($`C_{i+1}=C_{i+1}^{}`$) or rejected ($`C_{i+1}=C_i`$); the acceptance probability is given by

$$P_a=\mathrm{min}[1,\mathrm{exp}(\beta (E(C_{i+1}^{})E(C_i))],$$
(2)
in which $`\beta =1/(k_BT)`$ with Boltzmann constant $`k_B`$ and temperature $`T`$. The time scale is set such that in one unit of time, on average each spin is proposed to be flipped once. So in our system, in one unit of time we perform $`BL`$ Monte Carlo steps.

Our current interest in the Ising model with single-spin-flip dynamics stems from the fact that in finite systems, the configurations will occasionally switch between states in which the magnetization is either negative or positive, through an activated process. The dominant pathway at low temperatures consists of the formation of a single pair of closed interfaces in the shorter periodic direction (for $`BL`$), which perform a relative diffusive motion around the longer periodic direction and annihilate after meeting each other through the periodic boundary. A series of snapshots, illustrating this process, is presented in figure 1.

## III Theoretical framework

To study the behavior of times between magnetization reversals at temperatures below the critical one we may consider an ensemble of a large number of systems prepared in states with a magnetization close to the equilibrium value $`M_0`$ and study the rate at which these systems reach the value $`M=M_0`$, after which they are removed from the ensemble. Alternatively, the behavior of times between zero-crossings of the magnetization can be studied, as discussed below. The spin-flip dynamics described above may be represented by a master equation for the probability distribution $`P(\sigma )`$ of finding a system in the state $`\sigma `$ at time $`t`$. Due to the huge number of possible states this master equation cannot be solved analytically or even numerically for system sizes of practical interest. Therefore, as an approximation we assume that we may replace the exact master equation by an approximate master equation for the probability $`P(M,t)`$ of finding a system with magnetization $`M`$ at time $`t`$. The form of this equation is

$`{\displaystyle \frac{dP(M,t)}{dt}}`$ $`=`$ $`\mathrm{\Gamma }_{M,M+2}P(M+2,t)+\mathrm{\Gamma }_{M,M2}P(M2,t)`$ (3)
$``$ $`(\mathrm{\Gamma }_{M+2,M}+\mathrm{\Gamma }_{M2,M})P(M,t),`$ (4)
with $`\mathrm{\Gamma }_{M^{},M}`$ the transition rate from $`M`$ to $`M^{}`$. It ignores the fact that the actual transition rates will depend on the geometry of the state under consideration, with the idea that in typical cases spin states of very similar geometry will dominate the set of states with given $`M`$. Besides on $`M`$ and $`\beta `$ the transition rates will depend on the geometry of the system, that is on $`B`$ and $`L`$. To estimate the scaling behavior of these dependencies we make two further simplifying assumptions:

* First of all, requiring $`BL`$, we assume that states with $`M\pm M_0`$ consist of a single strip of opposite magnetization separated by two phase boundaries of length $`B`$ from the majority spin phase. The relevant changes of $`M`$ then will be caused by displacements of these boundaries due to flips of spins along them. The total number of spins available for this will be proportional to $`B`$, hence the transition rates should also be proportional to $`B`$ and independent of $`L`$.

This approximation will not be valid for values of $`M`$ that are too close to $`\pm M_0`$, as for these the opposite magnetization will typically be found in a closed cluster rather than in a strip. It will turn out though that the contributions from these $`M`$-values to the reversal frequency are very small for systems of reasonable length. Therefore using the approximation $`\mathrm{\Gamma }B`$ also here does not harm. Further we neglect the possibility of having more than one strip of opposite magnetization. When the temperature gets close to the critical one, when $`B`$ becomes small or when $`L`$ becomes very large, this may not be a good approximation. We will come back to this in section V. In our further theoretical treatment we will assume that we are in a situation where the single-strip approximation is justified.
* Further we neglect changes due to fluctuations of the magnetization caused by the growth and shrinkage of small clusters of the minority spin type. On average these do not contribute to magnetization reversals, so this approximation should be allowed. In case one considers the first passage frequency through zero magnetization there is a small effect, to which we will come back in section V too.

In order that the equilibrium distribution be a stationary solution of the master equation we impose the condition of detailed balance

$$\frac{\mathrm{\Gamma }_{M,M+2}}{\mathrm{\Gamma }_{M+2,M}}=\mathrm{exp}\left[\beta (F(M)F(M+2))\right],$$
(5)
where $`\beta F(M)=\mathrm{ln}P_{eq}(M)`$, with $`P_{eq}(M)`$ the equilibrium probability of magnetization $`M`$. In section IV B it will be explained how this condition may be combined with simulation results for interface diffusion to obtain approximations for all transition rates that lead to a good overall prediction of the reversal times.

The long-time reversal frequency as predicted by the master equation (4) follows as the largest eigenvalue of this equation, supplemented with an absorbing boundary at $`M=A`$, with either $`A=M_0`$ or $`A=0`$. The first choice corresponds to a real reversal of the magnetization, the second one to a first return to the value $`M=0`$. After this the system will have equal probabilities to actually reverse its magnetization or to return to the equilibrium magnetization value it came from, hence this return frequency should be twice the reversal frequency. The absorbing boundary condition is implemented by setting $`\mathrm{\Gamma }_{A2,A}`$ equal to zero.

The largest eigenvalue $`\nu `$ of $`\mathrm{\Gamma }_{M,M^{}}`$ in eq. (4), as well as the corresponding eigenvector $`P_0(M)`$, may be found by requiring that the net current away from magnetization $`M`$ assumes the value $`\nu P_0(M)`$. Using conservation of probability one easily checks that this may be expressed as

$`\mathrm{\Gamma }_{M+2,M}P_0(M)`$ $``$ $`\mathrm{\Gamma }_{M,M+2}P_0(M+2)=\nu {\displaystyle \underset{mM}{}}P_0(m),`$ (6)
$`\nu `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_{A,A2}P_0(A2)}{_{mA2}P_0(m)}}.`$ (7)
In eq. (6) $`P_0(m)`$ on the right hand side may be approximated, up to a normalization factor, by $`\mathrm{exp}(\beta F(m))`$, because the sum is dominated by the terms with small $`m`$-values, for which this approximation is excellent. This one may check in hindsight against the solution obtained. With this approximation the equation may be solved recursively for $`P_0(M)`$ in terms of $`P_0(A2)`$ for $`M=A4,A6,\mathrm{}`$, with the result

$`P_0(M)`$ $`=`$ $`{\displaystyle \underset{MmA2}{}}{\displaystyle \frac{\mathrm{\Gamma }_{A,A2}}{\mathrm{\Gamma }_{m+2,m}}}\mathrm{exp}\left[\beta (F(m)F(M))\right]`$ (8)
$`\times `$ $`{\displaystyle \frac{_{nm}\mathrm{exp}\left[\beta F(n)\right]}{_{n^{}A2}\mathrm{exp}\left[\beta F(n^{})\right]}}P_0(A2).`$ (9)
The summations over $`m`$ and $`M`$ are dominated by values of $`m`$ close to zero, combined with values for $`M`$ close to $`M_0`$, for which the sum over $`n`$ is basically independent of $`m`$. Therefore the reversal frequency may be obtained as

$`\nu `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_{A,A2}P_0(A2)}{_{MA2}P_0(M)}}`$ (10)
$`=`$ $`\left({\displaystyle \underset{m=M_0}{\overset{A2}{}}}{\displaystyle \frac{\mathrm{exp}\left[\beta F(m)\right]}{\mathrm{\Gamma }_{m+2,m}}}{\displaystyle \underset{nA2}{}}\mathrm{exp}\left[\beta F(n)\right]\right)^1.`$ (11)
The restriction of the summation over $`m`$ to values $`m>M_0`$ is needed to avoid large spurious contributions from $`m<M_0`$. The result in eq. (11) is well-known. It is usually derived by considering a state with a stationary current in which mass is inserted at a constant rate on one side (e.g. at $`M=B\times L`$ in our case) and taken out as soon as it reaches the absorbing boundary (see e.g. , section IV E). In that case the replacement of the sum over $`n`$ by a constant is exact.

## IV Simulations and results

In order to apply the above theoretical framework to magnetization reversal times in the Ising model, the two ingredients required are: (i) the equilibrium probability $`P_{eq}(M)`$ to find the system in a state with magnetization $`M`$; and (ii) the transition rates $`\mathrm{\Gamma }_{M^{},M}`$ from magnetization $`M`$ to $`M^{}`$. We obtain these two ingredients via two different computational approaches.

### A free energy landscape

For various values of $`L`$ and $`B`$ and various temperatures we make histograms of the distribution of the magnetization $`M=_i\sigma _i`$ and the energy. Since the probability that a state with energy $`E`$ occurs at an inverse temperature $`\beta `$ is proportional to the Boltzmann weight $`\mathrm{exp}(\beta E)`$, a histogram made at a certain temperature provides information about the probability distribution at nearby temperatures as well. We use the multiple histogram method  to combine the information from simulations at various temperatures. In this way, histograms for the magnetization over a wide range of temperatures can be obtained. As defined above, the free energy $`F(M)`$ is related to the probability that a certain magnetization $`M`$ occurs by $`P_{eq}(M)=e^{\beta F(M)}`$. The method of determining free energies of configurations constrained to some value of a coordinate along a pathway through phase space has been used extensively, for instance by Auer and Frenkel . In figure 2 the free energy as a function of the magnetization per spin $`m=M/N`$ is plotted at two different temperatures and three different values of $`L/B`$.

In order to compare these results with theory we note that the probability $`P_{eq}(0)`$ of finding the system in a state with $`M=0`$ can be obtained from eq. (4.23) of chapter V of  as follows. First note that this equation gives the following approximation for the partition function of a periodic anti-ferromagnetic system with one phase boundary , which is equivalent to a ferromagnetic system with an anti-periodic boundary:

$$Z_1=\frac{Z_0}{2}\underset{k=1}{\overset{2L}{}}\left(\frac{A(k)\sqrt{A^2(k)c^2}}{A(k)+\sqrt{A^2(k)c^2}}\right)^{B/2},$$
(12)
where

$`A(k)`$ $`=`$ $`(1+z^2)^22z(1z^2)\mathrm{cos}{\displaystyle \frac{\pi k}{L}},`$
$`c`$ $`=`$ $`2z(1z^2),`$
$`z`$ $`=`$ $`\mathrm{tanh}\beta J,`$
and $`Z_0`$ the partition function for a homogeneous system (cf. (4.20) of ). Neglecting the possible interaction between two phase boundaries (or interfaces) we conclude that the partition function $`Z_2`$ for a system with two interfaces is given by $`Z_2/Z_0=(1/2)(Z_1/Z_0)^2`$. The factor $`1/2`$ arises because otherwise we would count each configuration twice: interchanging the locations of the two interfaces does not give a different configuration. If we leave the positions of the two interfaces completely free the magnetization is distributed almost uniformly over all possible even values between $`M_0`$ and $`M_0`$. Therefore, fixing $`M`$ to 0 leads to a reduction of the partition function by a factor of $`M_0`$. Thus we arrive at the following result for $`P_{eq}(0)`$:

$`P_{eq}(0)`$ $`=`$ $`{\displaystyle \frac{Z_2}{M_0Z_0}}`$ (13)
$`=`$ $`{\displaystyle \frac{1}{8M_0}}\left[{\displaystyle \underset{k=1}{\overset{2L}{}}}\left({\displaystyle \frac{A(k)\sqrt{A^2(k)c^2}}{A(k)+\sqrt{A^2(k)c^2}}}\right)^{B/2}\right]^2.`$ (14)
For square systems this result has to be modified, because there are two ways to make a strip with opposite magnetization: the interfaces can lie in the horizontal as well as in the vertical direction. This gives an additional factor 2 in the equation for $`P_{eq}(0)`$.

### B interface diffusion coefficient

The second ingredient for the theoretical framework in section III consists of the transition rates $`\mathrm{\Gamma }_{M^{},M}`$ from magnetization $`M`$ to $`M^{}`$. As discussed there, we are actually mostly interested in the contribution to the transition rates that arises from the diffusion parallel to the longer periodic direction of interfaces that span the shorter one. To estimate this diffusion coefficient numerically, we study systems with anti-periodic boundary conditions: the spins in two neighboring sites $`i<N1`$ and $`j=i+1`$, or in $`i<NB`$ and $`j=i+B`$ are aligned if $`\sigma _i=\sigma _j`$, whereas the spins in two neighboring sites $`i=N1`$ and $`j=0`$, or in $`iNB`$ and $`j=i+BN`$, are aligned if $`\sigma _i=\sigma _j`$.

The interface location $`x(0)`$ is initially defined as the magnetization $`M(0)`$. As long as $`M[0.8N,0.8N]`$, steps $`\mathrm{\Delta }x`$ in the interface location are equal to changes $`\mathrm{\Delta }M`$ in the magnetization, i.e., $`x(t+\mathrm{\Delta }t)x(t)=M(t+\mathrm{\Delta }t)M(t)`$. Once the interface gets close to the anti-periodic boundary, the magnetization does not uniquely determine the interface location. As soon as $`M`$ is no longer in the interval $`[0.8N,0.8N]`$, we shift the anti-periodic boundary away from the interface over half the system size (which will bring the magnetization back in this range) and then continue.

In practice, we achieve this by switching from monitoring $`M`$ to monitoring $`M^{}=_{i=0}^{N/21}\sigma _i_{i=N/2}^{N1}\sigma _i`$, and $`\left|x(t+\mathrm{\Delta }t)x(t)\right|=\left|M^{}(t+\mathrm{\Delta }t)M^{}(t)\right|`$. The sign of the steps in $`x(t)`$ depends on the location of the interface: in the lower half this sign will be unchanged, whereas in the upper half it will be reversed. As soon as $`M^{}`$ leaves the interval $`[0.8N,0.8N]`$ we switch back to measuring steps in the original magnetization $`M`$. Also here, the sign of the steps in $`x(t)`$ depends on the location of the interface. Note that $`x`$ is not confined to the interval $`[N,N]`$. The diffusion coefficient $`D`$ is then obtained from the time-dependent interface location $`x(t)`$ as

$$D=\underset{t\mathrm{}}{lim}\left[\frac{(x(t)x(0))^2}{2t}\right].$$
(15)
We expect that the interface diffusion is independent of $`L`$, and that it grows linearly with $`B`$ since the number of sites where a spin flip moves the interface grows linearly with its width; however, for small $`B`$ corrections arise due to the periodicity (or helicity) of the boundaries. Consequently, we expect

$$D(B,L,\beta J)=g(\beta J)B+c.$$
(16)
For temperatures close to $`T_c`$ and small $`B`$, the system might occasionally contain more than a single interface. In that case, equation (15) would overestimate the diffusion coefficient of an interface. Since the additional free energy cost of interfaces increases linearly with $`B`$, this unwanted contribution to $`D`$ decreases exponentially with $`B`$. Taking this into account, our approach is to measure for various values of $`\beta J`$, $`B`$, and $`L`$ the diffusion coefficient $`D`$ via eq. (15). Next we determine the function $`g(\beta J)`$ in equation (16); the results are plotted in figure 3. Indeed we find that for increasing $`B`$ and at temperatures not too close to $`T_c`$ the diffusion coefficient rapidly becomes independent of $`L`$ within the ranges of $`L`$-values we considered. In our theoretical framework we then use as approximation for the jump rates

$$\mathrm{min}[\mathrm{\Gamma }_{M,M+2},\mathrm{\Gamma }_{M+2,M}]=\frac{g(\beta J)B}{2}.$$
(17)
This may be understood from the observations that on the one hand one has $`D=4\mathrm{\Gamma }`$, because the jumps in magnetization go by units of 2, on the other hand $`D`$ satisfies (16) with $`B`$ replaced by $`2B`$ because there are two interfaces. The above equation, in combination with eq. (5) and the free energy as a function of magnetization, specifies all the transition rates. Using (11) we can now predict the magnetization reversal times. The results are shown in Table I.

### C magnetization reversal times

We measure the magnetization every $`N`$ attempted spin flips. We look for events where the magnetization crosses from positive (more than half the spins up) to negative or vice versa between two consecutive measurements.

We make histograms of the times between two occurrences of magnetization reversal events. The histogram obtained in a system containing $`16\times 64`$ spins, at the temperature corresponding to $`\beta J=0.45`$, is shown in figure 4. For comparison we also show a histogram of the times measured between the first time the system reaches a free energy minimum, and the first time after this it reaches the other minimum. We will come back to this in our discussion.

The figure shows that at long times the decay function $`f(t)`$ behaves as $`f(t)\mathrm{exp}(t/\tau )`$. Here, we focus on the long-time behavior, and are specifically interested in the escape time $`\tau `$. We obtain this quantity via a fitting procedure, in which we ignore the data up to a time $`t_0`$, chosen such that $`f(t)`$ shows exponential time behavior for $`t>t_0`$. Then we determine the time $`t^{}`$ at which half of the remaining events have taken place. The escape time $`\tau `$ is then obtained from $`\tau =(t^{}t_0)/\mathrm{ln}(2)`$. Instead we could have made a linear fit of all the data points beyond $`t_0`$, but this makes no significant difference.

The resulting reversal times, for several system sizes and inverse temperatures, are presented in Table II and plotted in figure 5.

One sees by comparing the results of Tables I and II that in most cases the two agree within 20% for temperatures not too close to the critical temperature.

A commonly used first approximation to the description of the time scales of activated processes is Arrhenius’ law, which states that the typical time scale $`\tau `$ increases exponentially with the height of the (free) energy barrier $`\mathrm{\Delta }FF(M=0)F(M=M_0)`$, with a prefactor $`f`$ that depends in a mild way on temperature and system size,

$$\tau =f(\beta ,B,L)\mathrm{exp}(\beta \mathrm{\Delta }F).$$
(18)
If we straightforwardly use the free energy $`F(M)`$ as defined before we can obtain the heights of the free-energy barriers directly from the histograms of the magnetization distribution. To check the accuracy of this, we have plotted in figure 6 the ratio of $`\tau `$ and $`\mathrm{exp}(\beta \mathrm{\Delta }F)`$ as a function of inverse temperature $`\beta J`$, for various system sizes. Clearly, if the prefactor $`f`$ is assumed to be constant, this simple approximation fails to even predict the magnetization reversal times within an order of magnitude.

A little thought reveals that, in the present case, this way of determining the free energy barrier is not entirely satisfactory. The reason is that for magnetization values around $`M=M_0`$ the free energy increases with system size because the range of values through which $`M`$ typically fluctuates, increases as the square root of system size. As a consequence the probability of finding any specific $`M`$-value decreases. Around $`M=0`$ on the other hand, no similar effect occurs, due to the wide plateau of the free energy as function of $`M`$. As a result of this the decrease of probability of finding a given magnetization value for given relative position of the two interfaces is precisely compensated by the probabilities of finding this same magnetization value at different relative positions. Therefore it makes more sense in this case to define the Arrhenius factor in an alternative way as $`\mathrm{exp}(\beta F(M=0))`$. With this definition, where again $`\beta F(M)=\mathrm{ln}(P_{eq}(M))`$, it may be interpreted indeed as the equilibrium probability of finding two interfaces dividing the system into equal areas. But please notice that if, instead of a plateau, the free energy has a maximum of a width small compared to the range of magnetization fluctuations, the first definition of the Arrhenius factor is to be preferred.

With the second definition the expression for the decay time in terms of the Arrhenius factor becomes

$$\tau =f^{}(\beta ,B,L)\mathrm{exp}(\beta F(M=0)),$$
(19)
in which again one hopes that the prefactor $`f^{}`$ depends on temperature and system dimensions only in a mild way. On the basis of (11) and (17) we may conclude that for large systems $`f^{}`$ should simply be proportional to $`L`$ and independent of $`B`$.

Figure 7 shows $`f^{}`$ as a function of inverse temperature for several system sizes. For $`L=32`$ and 64 both the approximate independence of $`B`$ and the proportionality to $`L`$ are quite well-confirmed, especially if one takes into account that the plateau widths for these system sizes are notably less than $`L`$ (see figure 2). For $`L=16`$ the plateau becomes so narrow that the above predictions do not apply.

### D most slowly decaying mode

To obtain an estimate for the magnetization reversal times, we first estimated the most slowly decaying mode in eq. (9). Besides comparing the magnetization reversal times, we can also obtain the most slowly decaying eigenmode $`P_0(M)`$ directly from our Monte Carlo simulations by measuring the probability difference between the occurrence of a certain magnetization $`M`$ and the occurrence of the opposite magnetization $`M`$, averaged over a time scale comparable to the magnetization reversal time $`\tau `$. Since $`P_0`$ decays much slower than the other antisymmetric (in $`M`$) modes, it will give the dominant contribution to this probability difference. Figure 8 compares the most slowly decaying eigenmode as obtained with eq. (9) with the direct Monte Carlo measurements, for the $`16\times 64`$ system at $`\beta J=0.5`$.

## V Discussion

Our simulations confirm the global picture we sketched for the process of magnetization reversals in the Ising model with stochastic dynamics: a large cluster of opposite magnetization, originating through a fluctuation develops into a pair of interfaces. These interfaces diffuse around the system and annihilate, leaving the system in the oppositely magnetized phase. Quantitatively this process may be described to a good approximation as a diffusion process in a one-dimensional space with a coordinate describing the total magnetization in the system.

A few remarks should be made here.

* First of all it may seem remarkable that in almost all cases our theoretical prediction gives a shorter reversal time than the simulations. Since the theory neglects processes that enhance the reversal frequency, such as the formation of more than two simultaneous interfaces and spurious passages of $`M=0`$ (see below), one might expect an overestimate of the reversal time rather than an underestimate. The explanation of the latter comes primarily from that part of the process in which a growing cluster of opposite magnetisation transforms into a band around the cylinder. As the cluster will at first be typically circular in shape (if the temperature is not too low) it has to deform into more elliptical shape, with a longer interface than the band, before the latter can be formed. In figure 2 the signs of this are very clear in the form of little shoulders on the sides of the plateaus in the free energies. The transition rates $`\mathrm{\Gamma }_{M^{},M}`$ for the formation of these elliptical shapes will not be proportional to $`B`$ but much smaller, since this growth mainly occurs along the short sides of the cluster. The contributions from these terms in the denominator of (11) are therefore larger than estimated. Numerically they are important as long as the plateau in the free energy is not too wide. This then will give rise to a notable decrease of the reversal frequencies from our estimated values. Correcting for this in the theoretical expression would require better estimates of $`\mathrm{\Gamma }_{M^{},M}`$ in the free energy shoulders, which seem fairly complicated to obtain. Also the growth rates in the magnetization regions beyond the shoulders, where the growing clusters are even smaller, are overestimated by (17), but these are weighted less as the free energies there are smaller.
* The above considerations clearly reveal the possibility that the coordinate parametrizing the ”reaction path” (here the magnetization) changes non-monotonically along this reaction path; one could imagine that the growing cluster, before turning into a strip, typically decreases in size for a while. In reality this does not seem to happen, but if this were the case, parts of the reaction path with different cluster shape but equal magnetization would be lumped together, and the diffusion process along the cluster size coordinate would no longer correspond to the actual path taken by the cluster. Obviously in such a case the reaction coordinate has to be redefined in such a way that the new coordinate is monotonic along the reaction path indeed. But in complicated situations it may not always be clear what is a proper choice for such a coordinate.
* As a consequence of magnetization fluctuations in the bulk, passages of $`M=0`$ will be registered typically already a short while before the area between the interfaces reaches the value $`LB/2`$ corresponding to an equal division of the system between areas of positive and negative magnetization. The reason is that the typical time scale for fluctuations of the bulk magnetization is much shorter than that for interface diffusion, so for each location of the interfaces the whole range of accessible magnetization values typically will be scanned. Now in most cases this just will give rise to a negligible shift of the time at which the situation of equal areas (this is the physically relevant criterion) is reached, but occasionally it may happen that $`M=0`$ occurs, but the system returns to the pure state it came from without ever reaching equal areas. To check the importance of these events we may compare the average first passage time from $`M_0`$ to $`M_0`$ to that from $`\pm M_0`$ to $`M=0`$. If the effect of magnetization fluctuations is negligible the ratio of these should be 2, otherwise it ought to be larger. Note that for the systems considered here the effect of magnetization fluctuations on the first passage time from $`M_0`$ to $`M_0`$ is much smaller than that on the first passage time from $`M_0`$ to $`M=0`$, because the probability of returning to $`M_0`$ without reaching the pure state of positive magnetization, once $`M=M_0`$ has been reached, is extremely small. We have measured this ratio for several system sizes and temperatures, excluding those cases where the occurrence of multiple interfaces is likely, and found that the mean value is $`1.90`$, with a standard deviation of $`0.15`$.
* It is of interest to investigate how the reversal times depend on the system size parameters. From (11) we may conclude that for long systems (but not so long that multiple interface pairs will occur frequently) the reversal time will become independent of the system length $`L`$. For this one should notice first of all that for $`M`$-values on the plateau $`\mathrm{exp}(\beta F(M))`$ is proportional to $`1/L`$, because the number of ways a pair of interfaces may be placed such that the average magnetization equals $`M`$, is proportional to $`L`$. On the other hand the summation over $`M`$, which is dominated by $`M`$-values on the plateau, gives rise to a factor close to $`L`$. Hence to first approximation the reversal frequency $`\nu `$ is independent of $`L`$. Remarkably this independence of $`L`$ in fact is observed even better by the numerical results than by the theoretical estimates.

In section IV C it was noted already that the product $`f^{}`$ of the reversal frequency and the Arrhenius factor $`\mathrm{exp}(\beta F(0))`$ depends on $`L`$ and $`\beta J`$, but is approximately independent of the system width $`B`$. Notice that square systems may be included in this comparison without problem, as the extra factor of 2 due to the two possible orientations of the interfaces, is properly accounted for in the Arrhenius factor.
* We already indicated repeatedly that our theory may be applied only if the probability of having more than two interfaces around the cylinder may be neglected compared to the probability of having just two such interfaces. For this to be case one has to require $`Z_1/Z_01`$. For large enough systems this condition may be rewritten as $`L\mathrm{exp}(\beta \sigma B)`$, with $`\sigma `$ the surface tension of the interface. For large systems violation of this condition requires extreme aspect ratios, so under normal conditions it will be satisfied, unless the temperature is very close to the critical one, at which the surface tension vanishes.

We are presently applying the methods described here to a study of nucleation rates in metastable states. We hope to report on this before long.