format: text

text:
# Bosonic Excitations in Random Media

## I Introduction

Excitations in condensed matter systems with quenched disorder have been a subject of intense study during the last several decades. Historically, it has been fermionic excitations in random systems that have received most attention. The reason for this lies in part with the rapid development of experiments and theory involving mesoscopic conductors, where the effects of disorder in phase-coherent electron systems have been studied in great detail.Imry

It is, however, also of considerable interest to study random systems with bosonic excitations, and there is an extensive literature treating problems of this type, too. For instance, the propagation of phonons in glasses and of electromagnetic waves in media with random refractive index has long been a subject of active research,SJS ; PS  and trapping of light via scattering from disorder is a principle on which random lasers are based.Lasers  Other examples of bosonic excitations in random systems include vibrations of pinned elastic structures such as charge density waves,GL  magnons in diluted antiferromagnets and spin glasses,Harris ; WW ; HS ; AM ; Ginzburg ; Wan ; Chernyshev  and quasiparticles in superfluid liquid helium permeating a porous medium.Helium  To some extent, work on these problems has focussed on specific features of individual examples, and given less emphasis to a generic aspects than has been the case for disordered fermionic problems.

In this paper we emphasise just these generic aspects. In particular, we examine the relationship between universality classes identified for fermionic problems and models for bosonic excitations, as well as features that are specific to bosonic problems. In addition, we develop a framework for treating bosonic excitations which we use to calculate the density of states and localization properties of one dimensional, randomly pinned systems. Some of our results have been presented in short form elsewhere.GC

Our discussion is organised as follows. In Sec. II we review the symmetry classes established within random matrix theory for disordered fermionic systems. We recall in Sec. III the general form for a quadratic bosonic Hamiltonian and the Bogoliubov transformation required to diagonalise it, and show how it is useful to introduce an auxilliary problem with structure similar to that in the chiral symmetry class. In Sec. IV we discuss how various examples of systems with bosonic excitations fit into this general framework. Here we emphasise the distinction between excitations that are Goldstone modes and those that are not, and outline an established argument that leads in the latter case to the result $`\rho (\omega )\omega ^4`$. Then, as an illustrative example, we consider the random field $`XY`$ spin chain, using our formalism as the basis for a numerical study of the discrete system, presented in Sec. V, and giving an analytical treatment of the continuum limit in Sec. VI, recovering in both cases the behavior $`\rho (\omega )\omega ^4`$. In Sec. VII we consider the random field Heisenberg spin chain and some related but simpler models, which are of interest because disorder enters them in a more general way than for the $`XY`$ model. Finally, we summarise the relevant experimental situation in Sec. VIII and end with concluding remarks in Sec. IX

## II Symmetries and Disorder

To provide a context for our discussion of bosonic systems, we begin by setting out the symmetry classes that are recognized within random matrix theory for fermionic systems.Mehta ; Efetov  Models of non-interacting quasiparticles play an important role in the study of fermionic excitations in disordered conductors, insulators and superconductors. It has long been appreciated that the properties of these models are controlled by the discrete symmetries of the single particle Hamiltonians. Originally three symmetry classes for random Hamiltonians were identified.Wigner ; Dyson0  One class consists of random Hamiltonians that have time-reversal invariance but no Kramers degeneracy. In an appropriate basis the Hamiltonian $`H`$ is real, so that

$$H=H^{}.$$
(1)
Random Hamiltonians which obey Eq. (1) appear, for example, when studying disordered conductors without applied magnetic field. A second symmetry class consists of time-reversal invariant random Hamiltonians for particles with half-integer spin and hence Kramers degeneracy. In this case the time-reveral operation includes spin inversion, and invariance requires

$$H=\sigma _2H^{}\sigma _2.$$
(2)
(Here and in the following, $`\sigma _i`$ for $`i=1,2`$ or $`3`$ represent the conventional Pauli matrices, acting on a subspace identified by the context.) In order for Eq. (2) to be different in an essential way from Eq. (1), spin rotation invariance must be broken. Thus this case is of relevance for disordered conductors with spin-orbit coupling. A third symmetry class arises when the Hamiltonian has no discrete symmetries. Examples of these symmetry classes are provided by the three Wigner-Dyson random matrix ensembles.Wigner ; Dyson0

More recently, it has been recognized NS ; VZ ; Gade ; AZ  that there are seven additional classes of disordered fermionic Hamiltonians. These arise where there exists a special reference energy (taken to be zero in the following) for the system, and a symmetry operation which relates eigenstates in pairs. Three of these seven are referred to as chiral symmatery classes.NS ; VZ ; Gade  Hamiltonians for these classes can be put into the form

$$H=\left(\begin{array}{cc}0& Q\\ Q^{}& 0\end{array}\right)$$
(3)
where $`Q`$ is itself a matrix or an operator. They obey the symmetry condition

$$H=\sigma _3H\sigma _3.$$
(4)
This ensures that energy levels of such Hamiltonians appear in pairs $`\pm E`$, since if $`\psi `$ is an eigenfunction with energy $`E`$, then $`\sigma _3\psi `$ is an eigenfunction with energy $`E`$. The symmetry condition of Eq. (4), when combined with either Eq. (1), or Eq. (2), or neither, leads altogether to three chiral symmetry classes. Chiral Hamiltonians appear as tight-binding models with only off-diagonal disorderOE ; Eggarter ; Ziman ; Gade ; Brouwer  and in the problem of classical diffusion in a random medium.GB ; CW

The remaining four symmetry classes arise in the study of disordered superconductors with pairing treated in the mean field approximation. The Hamiltonian for such a problem has the Bogoliubov-de Gennes structure

$$H=\left(\begin{array}{cc}h& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& h^T\end{array}\right),$$
(5)
where the kinetic term $`h`$ is Hermitian, while the gap function $`\mathrm{\Delta }`$ is antisymmetric. This structure leads to another defining symmetry condition,

$$H=\sigma _1H^{}\sigma _1.$$
(6)
As for chiral Hamiltonians, the symmetry displayed in Eq. (6) ensures that energy levels of Bogoliubov-de Gennes Hamiltonians appear in pairs, $`\pm E`$. There are four symmetry classes of such Hamiltonians, according to whether or not the system has time-reversal and spin-rotation symmetry, making the total count ten.

A consequence of the conditions specified in Eqns. (4) and (6) is that statistical properties of energy levels and, in spatially extended systems, the associated eigenfunctions, are quite different in these additional symmetry classes near zero energy, compared to properties far from zero energy, or in the Wigner-Dyson symmetry classes.AZ

It is natural to ask whether this classification can be extended to problems involving non-interacting bosonic excitations or, equivalently, classical normal modes. At first sight, it might seem that quasiparticle statistics are unimportant in a non-interacting system. In one crucial respect, however, this is untrue, since stability of a system requires bosonic excitation energies to be positive, while for fermionic excitations stability is guaranteed by the Pauli exclusion principle, whatever the single-particle energy levels. This has two implications. First, energy zero emerges as a special point in the spectrum of bosonic systems, as it does for the additional fermionic symmetry classes discussed above. And second, one anticipates that matrix elements of bosonic Hamiltonians for random systems will have specific correlations, which ensure positivity of the spectrum. Thus, while the most general form for a quadratic bosonic Hamiltonian (see Eq. (7) below) is superficially similar to the Bogoliubov-de Gennes Hamiltonian, and while normal mode frequencies, like the eigenvalues of Eq. (5), appear in pairs $`\pm \omega `$, matrix elements must satisfy constraints in order that frequencies are real. Such a requirement is in stark contrast with the assumptions of statistical independence used in the construction of random matrix ensembles for fermionic systems. To stress the significance of this point, imagine a treatment of a disordered, interacting system which proceeds in two stages, by first finding the ground state and then calculating excitation energies within a harmonic approximation. The spirit of random matrix theory for fermionic systems is to divorce these two stages and approach the second one phenomenologically, choosing statistically independent matrix elements. By contrast, for bosonic systems it is clear that the two stages cannot be completely separated. Indeed, whereas for random matrix theory universal spectral properties follow largely from symmetry and are independent of the details of the matrix element distribution, we argue here that for bosonic excitations which are not Goldsone modes it is the requirements of stability and the ensuing correlations in the Hamiltonian for excitations that give rise to universal spectral properties.

## III Bosonic Hamiltonians

In this section we discuss the most general form for bosonic Hamiltonians and summarise the diagonalisation procedure, following standard lines (see, for example, Ref. Blaizot, ). We also emphasise the distinction between oscillator frequencies and stiffnesses. Finally, we show for oscillations about a stable ground state that it is natural to rewrite the Hamiltonian in terms of a chiral matrix.

### III.1 Stiffnesses and frequencies

The most general bosonic Hamiltonian can be written in the equivalent forms

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}\left[M_{ij}p_ip_j+K_{ij}q_iq_j+2C_{ij}q_ip_j\right]`$ (7)
$``$ $`{\displaystyle \frac{1}{2}}\left(\begin{array}{cc}𝐩& 𝐪\end{array}\right)\left(\begin{array}{cc}M& C\\ C^T& K\end{array}\right)\left(\begin{array}{c}𝐩\\ 𝐪\end{array}\right)`$
$``$ $`{\displaystyle \frac{1}{2}}\left(\begin{array}{cc}𝐚^{}& 𝐚\end{array}\right)\left(\begin{array}{cc}\mathrm{\Gamma }& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& \mathrm{\Gamma }^T\end{array}\right)\left(\begin{array}{c}𝐚\\ 𝐚^{}\end{array}\right)`$
Here, $`q_i`$ and $`p_i`$ are the coordinates and momenta of the oscillators, and $`a_i^{},a_i=\left(q_i\pm ip_i\right)/\sqrt{2}`$ are bosonic creation and annihilation operators. The matrices $`M`$ and $`K`$ are real and symmetric, while $`C`$ is an arbitrary real matrix. Equivalently, $`\mathrm{\Gamma }`$ is hermitian, while $`\mathrm{\Delta }`$ is symmetric. Physically, $`M`$ is the inverse mass matrix of the oscillators, $`K`$ is the matrix of spring constants, and couplings of the type represented by $`C`$ occur, for example, in spin systems. It is convenient to define the $`2N\times 2N`$ symmetric matrix

$$=\left(\begin{array}{cc}M& C\\ C^T& K\end{array}\right).$$
(8)
We are interested in frequencies of oscillators decribed by Eq. (7). From Hamilton’s equations of motion it follows that these frequencies are eigenvalues of the non-Hermitian matrix

$$^{}=\left(\begin{array}{cc}iC^T& iK\\ iM& iC\end{array}\right)\sigma _2,$$
(9)
They are real if the system is stable, in which case $`H`$ is bounded from below, the eigenvalues of $``$ are positive, and we can write $`=Q^TQ`$ with $`Q`$ real. In thse terms, to find frequencies we need to diagonalize the matrix $`\sigma _2Q^TQ`$, but its eigenvalues coincide with those of another matrix

$$\mathrm{\Omega }=Q\sigma _2Q^T,$$
(10)
which is explicitly Hermitian and, moreover, antisymmetric, so that frequencies come in pairs $`\pm \omega _i`$, with $`i=1\mathrm{}N`$.

If $`\mathrm{\Omega }`$ is interpreted as a random Hamiltonian, then according to the classification scheme discussed in Sec. II it belongs to one of the Bogoliubov-de Gennes classes (more precisely, to class D, see Ref. AZ, ). While this is indeed an indication that random oscillators behave in many ways like random fermionic Hamiltonians from one of the additional symmetry classes (having frequencies in pairs, with $`\omega =0`$ as a special point in the spectrum), the identification of $`\mathrm{\Omega }`$ as class D operator is not by itself necessarily helpful since $`\mathrm{\Omega }`$ does not have statistically independent matrix elements, but rather is constrained to have the form given in Eq. (10). Instead, we shall see that a link to matrices wih chiral symmetry proves more useful.

The computation of oscillator frequencies can equivalently be described as a Bogoliubov transformation for coordinates and momenta, specified by real matrices $`g`$ which obey

$$\sigma _2=g\sigma _2g^T$$
(11)
and tranform $``$ to $`gg^T`$. Diagonalizing $``$ using this transformation, the Hamiltonian $`H`$ of Eq. (7) takes the form

$$H=\underset{i}{}|\omega _i|a_i^{}a_i.$$
(12)
In addition to the frequencies, obtained as eigenvalues of $`\mathrm{\Omega }`$ or by Bogoliubov transformation of $``$, the eigenvalues of $``$ also have physical significance, for example as inverse static susceptibilities. We refer to them as stiffnesses, denoting them by $`\kappa _i`$, $`i=1,\mathrm{}2N`$. In general, there is no simple relationship between stiffnesses and frequencies, but several special cases provide important exceptions, as follows. Consider first $``$ with $`C=0`$ and $`M=𝟙`$, representing oscillations of particles, all with unit mass, connected by springs with spring constants $`K_{ij}`$. In this example, half of all the stiffnesses are equal to $`1`$, while the other half are the eigenvalues $`\kappa _i`$ of the matrix $`K`$; the frequencies and stiffnesses are related by

$$\omega _i=\pm \sqrt{\kappa _i}.$$
(13)
A second special case arises for magnon excitations in weakly disordered ferromagnets, which have a Hamiltonian of the form of Eq. (7) with $`M=K`$ and $`C=0`$. In that case, the stiffnesses are the eigenvalues of $`M`$ and $`K`$, and come in identical pairs. The frequencies are simply

$$\omega _i=\pm \kappa _i.$$
(14)
A final and important special case occurs when one stiffness, say $`\kappa _1`$, is much smaller than all others. In this regime, an approximate relation exists between the smallest frequency and the smallest stiffness:

$$\omega _1\pm \sqrt{\kappa _1}.$$
(15)
To derive this relation, a more detailed analysis of the structure of $``$ and $`\mathrm{\Omega }`$ is required. Since $``$ is a real symmetric matrix with positive eigenvalues, it can in general be represented as

$$=U\mathrm{\Lambda }^2U^T$$
(16)
where $`U`$ is an orthogonal matrix and $`\mathrm{\Lambda }_{ij}=\lambda _i\delta _{ij}`$ is the diagonal matrix whose eigenvalues are square roots of the eigenvalues of $``$. The corresponding matrix $`\mathrm{\Omega }`$ can be written as

$$\mathrm{\Omega }=\mathrm{\Lambda }U^T\sigma _2U\mathrm{\Lambda }.$$
(17)
Introducing an antisymmetric matrix $`A=U^T\sigma _2U`$, we have

$$\mathrm{\Omega }_{ij}=\lambda _iA_{ij}\lambda _j.$$
(18)
Suppose initially that one of the stiffnesses vanishes, so that $`\lambda _1=0`$. Then at least one frequency must vanish as well, since $`det\mathrm{\Omega }=det`$. In fact, because $`\mathrm{\Omega }`$ is antisymmetric and its frequencies come in opposite pairs, two frequencies are zero. The mathematical mechanism for this is clear from Eq. (18). First, since $`\mathrm{\Omega }_{1i}=\mathrm{\Omega }_{i1}=0`$ for all $`i`$, one of the eigenvalues of $`\mathrm{\Omega }`$ is $`\omega _1=0`$, with an eigenvector $`\psi _i^{(1)}=\delta _{1i}`$. Second, all other eigenvalues can be found by diagonalizing a smaller matrix, $`\mathrm{\Omega }_{ij}^{}`$ where $`2i,j2N`$. But $`\mathrm{\Omega }^{}`$ is an odd-dimensional antisymmetric matrix, and therefore has at least one zero eigenvalue, $`\omega _2`$, with an associated normalised eigenvector which we write as $`\psi _i^{(2)}`$, where $`\psi _1^{(2)}=0`$. Now treat small non-zero $`\lambda _1`$ using perturbation theory about this limit, with $`\lambda _1ϵ`$. The change in $`\mathrm{\Omega }`$ is

$`\delta \mathrm{\Omega }_{ij}=ϵ\left(\delta _{1i}A_{ij}\lambda _j+\lambda _iA_{ij}\delta _{1j}\right)(1\delta _{1i}\delta _{1j}).`$ (19)
Under this perturbation, the doubly degenerate eigenvalue $`\omega =0`$ of $`\mathrm{\Omega }`$ splits into $`\pm \omega _1`$, determined by diagonalizing the $`2\times 2`$ matrix

$$ϵ\left(\begin{array}{cc}0& _{i=2}^{2N}A_{1i}\lambda _i\psi _i^{(2)}\\ _{i=2}^{2N}A_{1i}\lambda _i\psi _i^{(2)}& 0\end{array}\right).$$
(20)
Hence

$$\omega _1=\pm \lambda _1\left|\underset{i=2}{\overset{2N}{}}A_{1i}\lambda _i\psi _i^{(2)}\right|\pm \sqrt{\kappa _1}$$
(21)
barring an accidental vanishing of the matrix element $`_{i=2}^{2N}A_{1i}\lambda _i\psi _i^{(2)}`$.

The usefulness of this result lies in the following. Consider a random system with bosonic excitations that are localized with a finite localization length at low frequency, and let the density of stiffnesses be $`d(\kappa )`$. We expect that each localization volume can be treated as an independent system, and that each low frequency excitation will have a frequency much smaller than that of other excitations in its own localization volume, and will therefore be associated with a single small stiffness. Applying Eq. (15), the density of excitation frequencies $`\rho (\omega )`$ for small $`\omega `$ is

$$\rho (\omega )=d(\omega ^2)\omega .$$
(22)
One interesting check of these conclusions is provided by random matrix theory. Consider an ensemble of real random matrices $`Q`$, with size $`N1`$ and probability distribution $`P\mathrm{exp}\left(Q^TQ\right)`$. Let $`\mathrm{\Delta }`$ be the typical magnitude of the eigenvalue of $`=Q^TQ`$ closest to zero. For $``$ generated in this way, one finds $`d(\kappa )\kappa ^{1/2}`$ at both $`\kappa \mathrm{\Delta }`$ and $`\kappa \mathrm{\Delta }`$. Computing the eigenvalue density of $`\mathrm{\Omega }=Q\sigma _2Q^T`$ using superintegrals, one finds thatMZ  $`\rho (\omega )`$ is independent of $`\omega `$ for $`\omega \mathrm{\Delta }`$, while $`\rho (\omega )\omega ^{1/3}`$ for $`\omega \mathrm{\Delta }`$. It is clear that in the regime $`\omega \mathrm{\Delta }`$, Eq. (22) is indeed applicable, while in the opposite regime it breaks down.

### III.2 Chiral Symmetry

To study stiffnesses of bosonic oscillators with $`=Q^TQ`$, it is advantageous to introduce $`\stackrel{~}{}`$, an auxiliary matrix in which $`Q`$ and $`Q^T`$ enter linearly.

$$\stackrel{~}{}=\left(\begin{array}{cc}0& Q\\ Q^T& 0\end{array}\right).$$
(23)
Obviously, the eigenvalues $`\lambda _i`$ of $`\stackrel{~}{}`$ are square roots of the stiffnesses $`\kappa _i`$ of $``$, and the matrix $`\stackrel{~}{}`$ plays the role of the square root of the original bosonic Hamiltonian $``$. On the other hand, the off-diagonal structure of $`\stackrel{~}{}`$ is the defining feature of the chiral symmetry class, discussed in Sec. II. While the direct implications of this connection are limited, because the elements of $`Q`$ are not independent random variables as would be the case in a random matrix ensemble, techniques originally developed for systems in this symmetry class will prove useful in our treatment of one-dimensional systems, as we describe in Sec. V and Sec. VI.

To study the frequencies of the oscillators, as opposed to their stiffnesses, a second auxiliary matrix

$$\stackrel{~}{}^{}=\left(\begin{array}{cc}0& Q\\ \sigma _2Q^T& 0\end{array}\right)$$
(24)
is helpful. We shall call matrices with this structure chiral bosonic matrices. The eigenvalue equation for $`\stackrel{~}{}^{}`$ can be written in the form

$$\left(\begin{array}{cc}0& Q\\ Q^T& 0\end{array}\right)\psi =\left(\begin{array}{cc}1& 0\\ 0& \sigma _2\end{array}\right)\sqrt{\omega }\psi ,$$
(25)
which of course inherits its structure from Hamilton’s equations.

The main difficulty in making use of these ideas is that, in practice, only $``$ is known initially and a method must be developed to find $`Q`$. Moreover, $`Q`$ is defined only up to a left multiplication by an arbitrary orthogonal matrix. For the introduction of $`Q`$ to be helpful, it will be important that it can be chosen to have a simple form, with, for example, only short range couplings. We shall show that this is indeed possible for a variety of problems.

## IV Bosonic Excitations: General Aspects

It is a feature of models for disordered fermionic systems from the chiral and Bogoliubov-de Gennes symmetry classes that their characteristic behavior appears only close to the reference energy, identified by the discrete symmetry of the Hamiltonian, while spectral properties at energies far from this are indistinguishable from those of the Wigner-Dyson classes. In a similar way, for bosonic excitations in disordered systems we expect it to be properties at low frequency that that are of particular interest. Examples of such excitations can be divided into two categories, according to whether or not they are Goldstone modes, associated with a broken continuous symmetry. In this section we summarize for each of these categories some of the previously established results. We also illustrate the introduction of a chiral Hamiltonian $`\stackrel{~}{}`$ for a one-dimensional phonon model.

### IV.1 Goldstone Modes

Important instances of Goldstone modes in disordered systems include acoustic phonons in glasses and alloys,SJS ; PS  and spin waves in isotropic, dilute ferromagnets and antiferromagnets, and isotropic spin glasses.Harris ; WW ; HS ; AM ; Ginzburg ; Wan ; Chernyshev  In disorder-free versions of these systems, low frequency excitations have long wavelength and are described by equations of motion that involve macroscopic properties of the system: density, elastic constants, magnetic susceptibility and spin stiffness. Disorder introduces local fluctuations in the values of these quantities, but one expects excitations to couple only to fluctuations averaged over a volume with linear dimensions set by a wavelength. Because of this, randomness only weakly affects low-frequency Goldstone modes, especially in higher dimensions for which the averaging is most effective. Such averaging is demonstrated by the low-frequency behavior of the excitation density $`\rho (\omega )`$, which above a critical dimension $`d_c`$ varies with the same power of $`\omega `$ as in the disorder-free system. It is also shown by localization properties of excitations: in one and two-dimensional systems, the localization length $`\xi (\omega )`$ diverges as $`\omega `$ approaches zero, while in higher-dimensional systems, all low-frequency states are extended. In this subsection we review the behavior of $`\rho (\omega )`$ and $`\xi (\omega )`$ for phonons in alloys and for spin waves in diluted antiferromagnets and spin glasses, considering the effect of weak disorder included in the relevant equation of motion. We also use the Hamiltonian for phonons in a one-dimensional disordered system as an illustration of the general mapping to chiral models, introduced in Sec. III.2.

Consider first a scalar version of a model for acoustic phonons, in which a mode with frequency $`\omega `$ has coordinate $`q(𝐫)`$ satisfying

$$\omega ^2q(𝐫)=c^2^2q(𝐫).$$
(26)
Suppose that the speed of sound, $`c`$ has random fluctuations in space, about an average value $`c_0`$, with only short-range correlations. Work on this and related problems is reviewed in Ref. PS, , and an early treatment of localization in this context was given in Ref. SJS, . The essentials for our purposes are as follows. First, the fluctuations in $`c`$, averaged over a $`d`$-dimensional volume of size set by the wavelength in the disorder-free system, decrease compared to $`c_0`$ with decreasing frequency as $`\omega ^{d/2}`$. Thus in this case the critical dimension is $`d_c=0`$, and for any $`d>d_c`$ the excitation density approaches the form found without disorder at low frequency. Moreover, from a calculation using the Born approximation, the Rayleigh scattering rate $`\tau ^1`$ varies as $`\tau ^1\omega ^{d+1}`$. (This dependence combines a factor of $`\omega ^2`$ from the frequency-dependence of coupling to disorder, and a factor of $`\omega ^{d1}`$ from the density of final states for scattering processes). In one dimension, disorder localizes with a localization length proportional to the mean free path, which here is $`c_0\tau `$, and soSJS

$$\xi (\omega )\omega ^2.$$
(27)
In two dimensions it is a familiar feature of electronic systems that the localization length is exponentially large in $`k_\mathrm{F}l`$ (where $`k_\mathrm{F}`$ is the Fermi wavevector and $`l`$ is the mean free path): the equivalent parameter for the phonon problem is $`\omega \tau `$, so thatSJS

$$\xi (\omega )\mathrm{exp}([\omega _0/\omega ]^2),$$
(28)
where $`\omega _0`$ is a disorder-strength dependent constant.

These results contrast interestingly with those for an antiferromagnet which has randomness generated either by site dilution (taken small enough that the system is above the percolation threshold) or by substitution of impurity spins with a magnetic moment different from that of the host spins.Harris ; Wan  In a discussion of the random antiferromagnet, it is useful to begin from the dispersion relation for spin waves in a two-sublattice ferrimagnet without disorder. At small wavevector $`k`$, this has the form

$$\omega ^2+(M_aM_b)\omega =c^2k^2,$$
(29)
where $`M_aM_b`$ is proportional to the difference between the two sublattice magnetic moments, which are taken to be oppositely aligned. Setting $`M_a=M_b`$, one recovers the usual dispersion relation in an antiferromagnet, with spin-wave speed $`c`$. For the random antiferromagnet, independent disorder on the two sublattices generates random fluctuations in the local value of $`M_aM_b`$ about a mean value of zero. Averaging these fluctuations over a volume of size set by the wavelength in the undiluted system gives a random variable with an amplitude that scales as $`\omega ^{d/2}`$. Because of this, disorder appears in Eq. (29) via the term $`(M_aM_b)\omega `$, which scales as $`\omega ^{d/2+1}`$. Above the critical dimension $`d_c=2`$, disorder is irrelevant in the sense that this term may be be neglected for small $`\omega `$ compared to $`\omega ^2`$. The value of the critical dimension is also apparent from a Born approximation calculation of the rate for scattering of spin waves by disorder,Wan  which yields $`\tau ^1\omega ^{d1}`$. Because of this, spin waves have a well-defined wavevector in the low frequency limit for $`d>d_c`$, since $`\omega \tau \mathrm{}`$ as $`\omega 0`$. Applying this approach below the critical dimension, for $`d=1`$ we see that Eq. (29) determines a relation between the lengthscale of an excitation, which we denote by $`1/k`$, and its frequency:

$$k^{d/2}\omega k^2$$
(30)
at small $`\omega `$. This implies

$$\rho (\omega )\omega ^{1/3}$$
(31)
and

$$\xi (\omega )\omega ^{2/3}.$$
(32)
We expect all these results for low-frequency behavior to be characteristic not only of disordered antiferromagnets, but also of spin glasses, since the two systems have in common the crucial feature of a magnetization density that is locally non-zero and random, but has zero average. A detailed treatment of excitations in spin glasses, however, is much more difficult than in weakly disordered antiferromagnets, because for spin glasses the ground state is generally disorder-dependent and unknown. Instead, the established approach is a hydrodynamic one,HS ; Ginzburg  which leads to linearly dispersing modes with a speed determined by the macroscopic spin-stiffness and susceptibility. In the light of our discussion, we expect these hydrodynamic results to be correct for $`d>d_c`$. By contrast, for one-dimensional systems, microscopic calculations are possible since frustration is absent and ground states can be determined. Results from computationsStinch  in one-dimensional models of $`\rho (\omega )`$ and $`\xi (\omega )`$ coincide with Eqns. (31) and (32), above. Similar calculations are also possible in higher dimensions for the Mattis model, which shares with one-dimensional models the features that frustration is absent and the ground state is known for each disorder configuration. In $`d=2`$ these yield $`\rho (\omega )\omega |\mathrm{log}(\omega )|`$, where the logarithm is characteristic of behavior at a critical dimension, and in $`d=3`$ they give $`\rho (\omega )\omega ^2`$, in agreement with hydrodynamic theory.Sher  A discussion of excitations in the Mattis model, building on the methods described in this paper, will be presented elsewhere.AG

Spin waves have also been investigatedBM  in the infinite-range Heisenberg spin glass: one of the main findings is a density of states that varies with frequency as $`\rho (\omega )\omega ^{\frac{3}{2}}`$. We have not been able to make contact between this result and the approaches described here.

As a next step, it is interesting to return to acoustic phonons in disordered systems and use the one-dimensional version of this problem to illustrate some of the methods set out in Sec. III. In fact, in this context the mapping to a chiral problem was exploited in celebrated early work by Dyson,Dyson  and also in calculations by Ziman.Ziman  Consider a one dimensional chain of particles with masses $`m_i`$, connected by nearest-neighbor springs with spring constants $`k_i`$, where $`m_i`$ and $`k_i`$ are random and positive. Let $`p_i`$ and $`q_i`$ be the momentum and displacement of the $`i`$-th particle The Hamiltonian is

$$H=\underset{i}{}\left[\frac{p_i^2}{2m_i}+\frac{k_i}{2}\left(q_iq_{i1}\right)^2\right],$$
(33)
It is convenient to make the canonical transformation $`p_i\sqrt{m_i}p_i`$ and $`q_iq_i/\sqrt{m_i}`$, giving

$$H=\underset{i}{}\left[\frac{p_i^2}{2}+\frac{k_i}{2}\left(\frac{q_i}{\sqrt{m_i}}\frac{q_{i1}}{\sqrt{m_{i1}}}\right)^2\right],$$
(34)
which is a particular case of the general bosonic Hamiltonian Eq. (7), in which $`M=𝟙`$ and $`C=0`$, so that eigenfrequencies are related to stiffnesses $`\kappa `$ that are eigenvalues of the matrix $`K`$, as in Eq. (13). Further discussion is easier in the continuum limit: replacing the index $`i`$ with a continuous coordinate $`x`$, the Hamiltonian becomes

$$H=𝑑x\left[\frac{p^2(x)}{2}+\frac{k(x)}{2}\left(\frac{d}{dx}\left[\frac{q(x)}{\sqrt{m(x)}}\right]\right)^2\right].$$
(35)
To find the stiffnesses, one must solve the eigenvalue equation $`Kq(x)=\kappa q(x)`$ with

$$Kq(x)=\frac{1}{\sqrt{m(x)}}\frac{d}{dx}k(x)\frac{d}{dx}\frac{q(x)}{\sqrt{m(x)}}.$$
(36)
Note that the operator $`K`$ must be positive definite, in order that the chain is stable. And indeed, defining $`a(x)=\sqrt{k(x)}`$ and $`b(x)=1/\sqrt{m(x)}`$, it may be expressed as a square, in the form $`K=Q^TQ`$ with $`Q=a(x)(d/dx)b(x)`$. Introducing a chiral Hamiltonian, as in Eq. (23), we then have

$$\left(\begin{array}{cc}0& a(x)\frac{d}{dx}b(x)\\ b(x)\frac{d}{dx}a(x)& 0\end{array}\right)\psi =\omega \psi .$$
(37)
At the equivalent point in his treatment, Dyson distinguishes between different possible choices for the form of disorder. In one case, termed type I, Eq. (37) is effectively replaced by Eq. (75) below (with $`V(x)=0`$), leading to the singularity of Eq. (61) in $`\rho (\omega )`$. An alternative, termed type II, retains instead Eq. (37) with only short-range correlations in $`a(x)`$ and $`b(x)`$, representing disorder that couples only weakly at small $`\kappa `$, because it is multiplied by spatial derivatives. This yieldsZiman  a constant $`\rho (\omega )`$ at small $`\omega `$, in agreement with the general arguments set out following Eq. (26). ZimanZiman  has given a detailed discussion in this context of the consequences of different types of disorder.Comment

### IV.2 Non-Goldstone Low Energy Excitations

Without Goldstone modes, the very existence of the low-lying excitations on which we have focussed our attention is not guaranteed. In fact, as seems first to have been appreciated in the context of atomic vibrations in glasses,IK  disorder itself may provide a route to a gapless spectrum. The essential ingredients are that the ground state should depend on the disorder realization, and that excitations at low frequency should be localized by disorder. Then it is reasonable to consider excitations within each localization volume separately, and to expect disorder configurations that support low frequency excitations to occur with finite probability. Roughly speaking, these excitations occur in regions where the ground state configuration is unusally sensitive to small changes in the disorder. In this section we summarise an established approach to this phenomenon of disorder-generated low frequency excitations, which concentrates on a single coordinate and its conjugate momentum. In subsequent sections we apply the formalism developed in this paper to study the phenomenon more generally.

Following Il’in and Karpov,IK  consider a one-dimensional anharmonic oscillator with momentum $`p`$, mass $`m`$, coordinate $`q`$ and potential $`U(q)`$. The Hamiltonian is

$$H=\frac{p^2}{2m}+U(q).$$
(38)
Choosing $`U(q)`$ to be a smooth random function, we wish to study the frequency of small-amplitude oscillations about the absolute minimum in $`U(q)`$ and, specifically, the probability distribution of this frequency. To give a more precise meaning to the notion of a smooth random function, we expand it in Taylor series

$$U(q)=\underset{n=1}{\overset{\mathrm{}}{}}a_n\frac{q^n}{n!},$$
(39)
where, to fix the zero of energy, we set $`U(0)=0`$. We take the $`a_n`$ to be random with joint probability distribution $`P(a_1,a_2,\mathrm{})`$. We shall assume that $`P(a_1,a_2,\mathrm{})`$ is free of zeros and divergences, but its detailed form will not be important.

Suppose that $`U(q)`$ has a minimum at $`q=q_0`$. In order to discuss excitations of the oscillator, we first Taylor expand $`U(q)`$ around $`q=q_0`$, writing

$$U(q)=\underset{n=2}{\overset{\mathrm{}}{}}b_n\left(\frac{(qq_0)^n}{n!}\frac{(q_0)^n}{n!}\right).$$
(40)
The probability distribution of the coefficients $`b_n`$ is related to that for the $`a_n`$ by

$$P(q_0,b_2,b_3,\mathrm{})=\left|det\frac{(a_1,a_2,a_3,\mathrm{})}{(q_0,b_2,b_3,\mathrm{})}\right|P(a_1,a_2,a_3,\mathrm{}).$$

Evaluating the Jacobian, we find

$$P(q_0,b_2,b_3,\mathrm{})=|b_2|P(a_1,a_2,a_3,\mathrm{}).$$
(41)
The probability distribution of $`b_2`$, the curvature of the potential $`U(q)`$ at a turning point, is hence

$$P(b_2)=|b_2|𝑑q_0𝑑b_3𝑑b_4\mathrm{}P(a_1,a_2,a_3,\mathrm{}).$$
(42)
Under the assumption that $`P(a_1,a_2,\mathrm{})`$ is free of zeros and divergences, this integral remains finite as $`b_20`$, and so for small $`b_2`$ we have $`P(b_2)b_2`$.

Small amplitude oscillations about $`q_0`$ have a frequency $`\omega \sqrt{b_2}`$, and it then follows that the probability distribution of these oscillation frequencies varies for small $`\omega `$ as

$$\rho (\omega )\omega ^3.$$
(43)
It is a further restriction to demand that a minimum at $`q_0`$ is the global minimum of $`U(q)`$. A full treatment of this constraint would be difficult but is fortunately not necessary: the crucial condition is that $`U(q)`$ should have no nearby minima deeper than the one at $`q_0`$. For that it is sufficient to truncate the expansion

$`U(q)=`$ $`U(q_0)+{\displaystyle \frac{b_2}{2}}(qq_0)^2+{\displaystyle \frac{b_3}{6}}(qq_0)^3`$ (44)
$`+{\displaystyle \frac{b_4}{24}}(qq_0)^4+\mathrm{}`$
at $`𝒪(qq_0)^4`$. After truncation, we require $`|b_3|<\sqrt{3b_2b_4}`$ in order that $`U(q_0)U(q)`$ for all $`q`$. This further suppresses the probability density for small curvatures, giving

$$P(b_2)=b_2𝑑q_0𝑑b_4\mathrm{}_{\sqrt{3b_2b_4}}^{\sqrt{3b_2b_4}}𝑑b_3P(a_1,a_2,\mathrm{})b_2^{\frac{3}{2}}.$$

In turn, this brings the frequency distribution at small $`\omega `$ to the form

$$\rho (\omega )\omega ^4.$$
(45)
One expects that the result of including higher order terms in Eq. (44), and of ensuring that no more distant minima are lower than the one at $`q_0`$, will be to change the constant of proportionality but not the power in Eq. (45).

Clearly, a serious limitation of this discussion it that it is limited to a system with a single coordinate and conjugate momentum. As with our closely-related discussion of random matrices, preceding Eq. (22), we expect the result to apply quite generally, provided excitations are localized with a finite localization length at low frequency. In these circumstances, the coordinate $`q`$ is interpreted as being the relevant degree of freedom for a low-frequency excitation within a localization volume. One of our objectives in Sec. V and Sec. VI is to provide detailed evidence for the more general relevance of Eq. (45).

## V Discrete 1D Random Field XY model

In this section and Sec.VI we study the excitations of the one-dimensional, classical $`XY`$ spin chain in a random field. This model provides a simple but non-trivial example of a system with bosonic excitations which are not Goldstone modes, and has been discussed previously in Refs. AR,  and Fogler, . In what follows we apply the general approach described in Sec. III. First, we set out definitions and write down the Hamiltonian $``$ for small amplitude excitations. Second, we find a local $`Q`$ which satisfies $`=Q^TQ`$. As a result, we map our problem onto a well-known one, involving a one-dimensional chiral Hamiltonian. We next review established results for such one-dimensional chiral Hamiltonians, which can be in one of two regimes, depending on the details of their disorder distribution. We use the mapping to obtain $`\rho (\omega )`$ both numerically, and (in Sec. VI) analytically.

### V.1 Definitions

The Hamiltonian for the random field $`XY`$ spin chain is

$$H=\frac{1}{2}\underset{i=1}{\overset{N}{}}\mathrm{\Pi }_i^2\underset{i}{}\mathrm{cos}(\varphi _i\varphi _{i+1})\underset{i}{}h_i\mathrm{cos}(\varphi _i\chi _i)$$
(46)
Here $`\mathrm{\Pi }_i`$ are momenta conjugate to the spin angles $`\varphi _i`$, exchange energy is represented by $`\mathrm{cos}(\varphi _i\varphi _{i+1})`$, and $`h_i`$ and $`\chi _i`$ are the amplitude and phase of a random field. It is convenient to introduce the notation $`I(\varphi )\mathrm{cos}(\varphi )`$ and $`h_i(\varphi )h_i\mathrm{cos}(\varphi \chi _i)`$.

While the kinetic energy is quadratic in the momenta $`\mathrm{\Pi }_i`$, the potential energy is strongly anharmonic in the coordinates $`\varphi _i`$. We want to find the ground state spin configuration $`\varphi _i^0`$ and the frequencies of oscillations about that ground state. The ground state configuration satisfies

$$\frac{H}{\varphi _i}|_{\varphi =\varphi ^0}=0.$$
(47)
Expanding $`H`$ about $`\varphi _i^0`$ to quadratic order, Eq. (46) reduces to an expression of the general form given in Eq. (7), and specified by the matrices $`C`$, $`M`$ and $`K`$. In this case, $`C=0`$ and $`M=𝟙`$. The symmetric matrix of spring constants, $`K_{ij}=^2H/\varphi _i\varphi _j|_{\varphi =\varphi ^0}`$ is tridiagonal, with no-zero entries

$$K_{ii}=_\varphi ^2I(\varphi _i^0\varphi _{i1}^0)+_\varphi ^2I(\varphi _i^0\varphi _{i+1}^0)+_\varphi ^2h_i(\varphi _i)$$
(48)
and

$$K_{i,i+1}=K_{i+1,i}=_\varphi ^2I(\varphi _i^0\varphi _{i+1}^0).$$
(49)
As discussed in Sec. III, with $`M`$, $`C`$ and $`K`$ of this form, the excitation frequencies of the system are $`\omega _i=\pm \sqrt{\kappa _i}`$, where the stiffnesses $`\kappa _i`$ are the eigenvalues of $`K`$. Our tasks, then, are the linked ones of determining $`\varphi _i^0`$ and diagonalizing $`K`$.

### V.2 Mapping onto a chiral problem

Since $`K_{ij}`$ is a real positive matrix, it can be written as the square of another real matrix $`Q`$, in the form $`K=Q^TQ`$. Our strategy is to find $`Q`$ and then study the related chiral Hamiltonian

$$\stackrel{~}{}=\left(\begin{array}{cc}0& Q\\ Q^T& 0\end{array}\right).$$
(50)
As the frequencies coincide with the square roots of the eigenvalues of $`K`$, we will use $`\omega `$ to denote the eigenvalues of $`\stackrel{~}{}`$ in this and the following section.

Because $`K`$ is tridiagonal, the matrix $`Q`$ can be chosen bidiagonal, with non-zero elements $`Q_{ii}A_i`$ and $`Q_{i,i1}B_{i1}`$ which satisfy

$$K_{ii}=A_i^2+B_i^2,K_{i,i+1}=A_{i+1}B_i.$$
(51)
To show this and find $`A_i`$ and $`B_i`$, it is helpful to use the idea of a partial energy, first introduced in this context by Feigelman: Fei  $`_i(\varphi _i)`$ is defined to be the ground state energy of a subsystem which consists of sites $`ji`$, considered as a function of the orientation of the boundary spin, $`\varphi _i`$. Formally

$$_i(\varphi _i)=\underset{\varphi _j,j<i}{\mathrm{min}}\underset{j<i}{}\left[I(\varphi _j\varphi _{j+1})+h_j(\varphi _j)\right]+h_i(\varphi _i).$$
(52)
It satisfies the recursion relation

$$_i(\varphi )=\underset{\psi }{\mathrm{min}}\left[I\left(\varphi \psi \right)+_{i1}(\psi )\right]+h_i(\varphi ).$$
(53)
Now let the value of $`\psi `$ which results from the minimization in Eq. (53) be $`\psi _0(\varphi )`$, and define $`\varphi _{i1}`$ as a function of $`\varphi _i`$ by $`\varphi _{i1}(\varphi _i)\psi _0(\varphi _i)`$. With this notation, the condition that the right side of Eq. (53) is at a minimum takes the form

$$_\varphi I(\varphi _i\varphi _{i1})+_\varphi _{i1}(\varphi _{i1})=0.$$
(54)
By differentiating Eq. (54) with respect to $`\varphi _i`$, remembering that $`\varphi _{i1}`$ is a function of $`\varphi _i`$ in the sense described above, we find

$$\frac{d\varphi _{i1}}{d\varphi _i}=\frac{_\varphi ^2I(\varphi _i\varphi _{i1})}{_\varphi ^2I(\varphi _i\varphi _{i1})+^2_{i1}(\varphi _{i1})}.$$
(55)
In addition, we differentiate (53) twice with respect to $`\varphi _i`$, again remembering that $`\varphi _{i1}`$ is a function of $`\varphi _i`$, to find

$`_\varphi ^2h_i(\varphi _i)=_\varphi ^2_i(\varphi _i)`$ (56)
$`_\varphi ^2I(\varphi _i\varphi _{i1})`$ $`{\displaystyle \frac{_\varphi ^2_{i1}(\varphi _{i1})}{_\varphi ^2I(\varphi _i\varphi _{i1})+^2_{i1}(\varphi _{i1})}}.`$
This allows us to solve Eqns. (48), (49) and (51), obtaining

$`A_i^2={\displaystyle \frac{\left[_\varphi ^2I(\varphi _i^0\varphi _{i1}^0)\right]^2}{_\varphi ^2I(\varphi _i^0\varphi _{i1}^0)+_\varphi ^2_{i1}(\varphi _{i1}^0)}}`$ (57)
and

$`B_i^2=_\varphi ^2I(\varphi _{i+1}^0\varphi _i^0)+_\varphi ^2_i(\varphi _i^0),`$ (58)
which completes the derivation of $`Q`$.

Let us examine the chiral Hamiltonian $`\stackrel{~}{}`$, Eq. (50). By rearranging its rows and columns it can be put into the form

$$\stackrel{~}{}=\left(\begin{array}{cccccccc}0& B_1& 0& 0& 0& \mathrm{}& 0& 0\\ B_1& 0& A_2& 0& 0& \mathrm{}& 0& 0\\ 0& A_2& 0& B_2& 0& \mathrm{}& 0& 0\\ 0& 0& B_2& 0& A_3& \mathrm{}& 0& 0\\ 0& 0& 0& A_3& 0& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}& 0& B_N\\ 0& 0& 0& 0& 0& \mathrm{}& B_N& 0\end{array}\right)$$
(59)
which is familiar as a one-dimensional tight binding model with only off-diagonal disorder,Eggarter ; Brouwer  and is also referred to as a random chiral one-dimensional Hamiltonian.

### V.3 Established results for one-dimensional chiral problems and implications for $`XY`$ model

There has been extensive previous work on one dimensional models of the type represented by Eq (59), with disorder in the $`A_i`$ and $`B_i`$ which is uncorrelated and chosen to have a simple, known distribution. The results serve as a guide for our calculations, and we summarise them here. Parameterise the matrix elements $`A_i`$ and $`B_i`$ as

$$A_i=1+a+\delta A_i,B_i=1a+\delta B_i$$
(60)
and take $`\delta A_i`$ and $`\delta B_i`$ to be independent Gaussian random variables with zero mean and standard deviation $`\sigma `$. Behavior at small $`\omega `$ is very different according to whether or not $`a`$ is zero.

For $`a=0`$ and $`\omega \sigma `$, the density of states has a singularityEggarter  at $`\omega =0`$ of a type first obtained in a related problem by Dyson,

$$\rho (\omega )\frac{1}{\omega \left|\mathrm{log}^3\omega \right|},$$
(61)
and the localization length of these states diverges for $`\omega 0`$ as

$$\xi (\omega )=|\mathrm{log}(\omega )|.$$
(62)
By contrast, for $`a0`$ (sometimes referred to as the staggered regime Brouwer ) the density of states varies asOE

$$\rho (\omega )\omega ^\beta ,\omega a,$$
(63)
with a power $`\beta `$ that depends on $`a`$ and $`\sigma `$. The localization length is finite in the small $`\omega `$ limit, and varies for small $`a`$ as $`\xi |a|^1`$.

While the $`A_i`$ and $`B_i`$ which arise in our treatment of the $`XY`$ model have specific correlations not present in the chiral problems studied previously, some comparisons are useful. In particular, considering for simplicity weak disorder, $`h_i1`$, the system arising from the $`XY`$ model turns out to be in the staggered regime. To see this, note that for weak disorder, $`|\varphi _i^0\varphi _{i1}^0|1`$, so that $`_\varphi ^2I(\varphi _i^0\varphi _{i1}^0)1`$,

$$A_i1\frac{1}{2}_\varphi ^2_{i1}(\varphi _{i1}^0),B_i1+\frac{1}{2}_\varphi ^2_i(\varphi _i^0)$$
(64)
and hence $`a=_\varphi ^2/2`$. We argue in Sec. VI that $`_\varphi ^2>0`$. We therefore expect low frequency states to be localized with localization length $`\xi _\varphi ^2^1`$, and with a power law density, as in Eq. (63). Quite separately, if states are localized, the assumptions that led to result $`\rho (\omega )\omega ^4`$ in Sec. IV are justified, and so we expect the exponent $`\beta =4`$. We postpone further analytical work to Sec. VI, and first treat the problem numerically.

### V.4 Numerical study of the one-dimensional random field $`XY`$ model

Our numerical procedure involves several steps. First, for a system of length $`L`$, we construct the function $`_i(\varphi )`$ for each $`i`$ by iterating Eq. (53) from $`i=1`$ to $`i=L`$. Then we determine $`\varphi _i^0`$, iterating from $`i=L`$ to $`i=1`$ and using the fact that for each $`\varphi _i^0`$, $`\varphi _{i1}^0`$ minimises the right-hand side of Eq. (53). Knowing the ground state spin configuration, we compute the matrix elements appearing in the chiral Hamiltonian Eq. (59), using Eqns. (57) and (58). Finally, we employ the transfer matrix technique developed specifically for such Hamiltonians in Ref. Eggarter,  to find the integrated density of states. For the random field we choose a uniform distribution of $`[h_i\mathrm{cos}(\chi _i),h_i\mathrm{sin}(\chi _i)]`$ over a disc of radius $`D`$, independently for each $`i`$. We find that it is sufficient for each disorder strength $`D`$ to study a single realization in a system of length $`L=10^6`$.

Our results for the integrated density of stiffnesses, $`N(\omega )=_0^\omega 𝑑\omega ^{}\rho (\omega ^{})`$, are shown in Fig 1. Behavior at large $`\omega `$ approaches that in the disorder-free system, as indicated. At small $`\omega `$ we expect the power law

$$N(\omega )\omega ^{(\beta +1)}.$$
(65)
This form and also the value $`\beta =4`$ are both already apparent in Fig 1, as was also the case, at lower precision, in earlier numerical results of FoglerFogler  for much shorter chains. It is more instructive, however, to observe that for a power-law density, dimensional analysis fixes the dependence on disorder strength $`D`$ to be $`\rho (\omega )D^{2\beta /3}\omega ^\beta `$. In Fig. 2 we plot $`D^{8/3}N(\omega )`$ as a function of $`\omega `$, showing collapse of data at small $`\omega `$ for three different disorder strengths, and power-law behavior with $`\beta =4`$.

We conclude that, while the exponent $`\beta `$ appearing in the density of states is, for a generic chiral problem, disorder-dependent and hence non-universal, the particular disorder generated in the mapping from the random spin chain, following Eqns. (53), (57) and (58), has a specific distribution and the correlations necessary to produce a universal density of states, $`\rho (\omega )\omega ^4`$. In the next section we present the analytic derivation of this result.

## VI Analytic treatment of the random field $`XY`$ model in the continuum limit

In this section we continue our examination of excitations in the one-dimensional random field $`XY`$ model, using the continuum limit to make analytic progress. We find both the localization length as a function of disorder strength and the density of states as a function of frequency, confirming heuristic arguments and numerical results given above.

### VI.1 Mapping to a chiral problem

The continuum limit of the random field $`XY`$ spin chain is reached at weak disorder, $`h_i1`$. In this limit it is possible to replace the discrete index $`i`$ with a continuous variable $`x`$, and the Hamiltonian of Eq. (46) becomes

$$H=𝑑x\left[\frac{1}{2}\mathrm{\Pi }^2+\frac{1}{2}\left(_x\varphi \right)+h(\varphi (x),x)\right],$$
(66)
with $`h(\varphi (x),x)=h(x)\mathrm{cos}\left(\varphi (x)\chi (x)\right)`$. As in the discrete case, we are interested in the configuration of spins $`\varphi _0(x)`$ that minimizes the potential energy. This configuration satisfies

$$_x^2\varphi _0(x)+_\varphi h(\varphi _0(x),x)=0$$
(67)
and the amplitudes of normal mode excitations about the ground-state $`\varphi _0`$ obey

$$\left[_x^2+_\varphi ^2h(\varphi _0(x),x)\right]\psi =\omega ^2\psi .$$
(68)
Just as in the discrete case, because this equation describes deviations from a minimum, all normal modes have positive stiffness $`\kappa =\omega ^2`$. As a result, the operator appearing in Eq. (68) can be written as a square. We set

$$\left[\frac{d^2}{dx^2}+_\varphi ^2h(\varphi _0(x),x)\right]=Q^TQ,$$
(69)
where $`Q=d/dx+V(x)`$. In other words, we require a function $`V(x)`$, which we term the chiral potential, that satisfies

$$\frac{dV(x)}{dx}+V^2(x)=_\varphi ^2h(\varphi _0(x),x).$$
(70)
In order to understand properties of the chiral potential, it is useful to introduce $`(\varphi ,x)`$, a continuum version of the partial energy,Fei  defined following Eq. (52) by

$$(\phi ,x)=\underset{\varphi (x)=\phi }{\mathrm{min}}_0^x𝑑y\left[(_y\varphi )^2+h(\varphi (y),y)\right].$$
(71)
Similarly, the continuum version of Eq. (53) is

$$_x+\frac{1}{2}\left(_\varphi \right)^2=h(\varphi ,x),$$
(72)
which can be thought of as a Hamilton-Jacobi equation for the action $``$ of a particle with coordinate $`\varphi `$ moving as a function of time $`x`$ in the time-dependent potential $`h(\varphi (x),x)`$. In addition, the continuum version of Eq. (54) relates the ground state configuration to this action, via

$$\frac{d\varphi _0(x)}{dx}=_\varphi (\varphi _0(x),x).$$
(73)
It is easy to check that by writing

$$V(x)_\varphi ^2(\varphi _0(x),x)$$
(74)
we solve Eq. (70). Thus the chiral potential $`V(x)`$ may be expressed simply in terms of the dependence of the ground state energy of a half-system on the boundary spin, $`\varphi (x)`$: the half-system has coordinate $`y`$ taking values in the range $`0yx`$. With $`V(x)`$ in hand, we wish to study the continuum version of Eq. (59): the chiral operator

$$\stackrel{~}{}=\left(\begin{array}{cc}0& Q\\ Q^T& 0\end{array}\right)=\left(\begin{array}{cc}0& \frac{d}{dx}+V(x)\\ \frac{d}{dx}+V(x)& 0\end{array}\right),$$
(75)
which has eigenvalues $`\pm \omega `$ that are the square roots of those appearing in Eq. (68).

### VI.2 Treatment of one-dimensional chiral problems

As with the lattice version, discussed in Sec.V.3, the spectral properties of $`\stackrel{~}{}`$, Eq. (75), have been studied extensively with simple choices for the probability distribution of $`V(x)`$. Behavior is as summarised for the lattice version in Sec.V.3; a particularly detailed study is given in Ref. Comtet, . Here, for completeness we sketch the derivation of the result that is of most importance for our work: the density of states at low frequency in the staggered regime, where $`V(x)>0`$ for the continuum system plays the same role as $`a>0`$ for the lattice model.

Following Ref. Comtet, , consider the coupled first-order differential equations

$$\left(\begin{array}{cc}0& \frac{d}{dx}+V(x)\\ \frac{d}{dx}+V(x)& 0\end{array}\right)\left(\begin{array}{c}\psi _1(x)\\ \psi _2(x)\end{array}\right)=\omega \left(\begin{array}{c}\psi _1(x)\\ \psi _2(x)\end{array}\right)$$
(76)
for $`x>0`$, with boundary conditions $`\psi _1(0)=0`$, $`\psi _2(0)=1`$. From the node counting theorem, the integrated density of states is given by the density of zeros of $`\psi _1(x)`$ per unit length. Introducing the parametrization $`\psi _1(x)=\rho (x)\mathrm{sin}\theta (x),\psi _2(x)=\rho (x)\mathrm{cos}\theta (x)`$, one has

$$\frac{d\theta (x)}{dx}=\omega V(x)\mathrm{sin}\left(2\theta (x)\right).$$
(77)
Thus we require the average rate of increase in phase with length, $`d\theta (x)/dx`$. To find this at small $`\omega `$, note that, by assumption, $`V(x)`$ is mainly positive, and for positive $`V(x)`$ Eq. (77) has stable fixed points close to $`\theta =n\pi `$, with $`n`$ integer. Rare fluctuations of $`V(x)`$ which are negative for a long interval in $`x`$ allow $`\theta (x)`$ to grow, evolving with increasing $`x`$ from one such fixed point to the next. Suppose $`\theta (x)`$ leaves the vicinity of one fixed point at $`x=x_1`$ and arrives in the vicinity of the next at $`x=x_2`$. For $`x_1<x<x_2`$ we neglect $`\omega `$ in Eq. (77) and obtain

$$\frac{d\theta }{\mathrm{sin}\left(2\theta \right)}=𝑑xV(x).$$
(78)
We estimate the integral on the left hand side by noting that the most important contribution comes from the end points, where $`\theta (x_1)n\pi `$, $`\theta (x_2)(n+1)\pi `$ and $`\mathrm{sin}\left(2\theta \right)\omega `$. For $`\theta `$ to increase by $`\pi `$ we therefore require a fluctuation in $`V(x)`$ which is sufficiently negative and extends over a sufficiently large interval in $`x`$ that

$$\mathrm{log}\left(\omega \right)>_{x_1}^{x_2}𝑑xV(x).$$
(79)
Let $`P(\omega )`$ be the probability per unit length for such a fluctuation to occur: the integrated density of states, $`_0^\omega \rho (\omega )𝑑\omega `$, is simply $`P(\omega )`$. It is natural to expect this probability to be exponentially small in $`\mathrm{log}(\omega )`$ for small $`\omega `$, so that, introducing a constant $`\alpha `$,

$$P(\omega )\mathrm{exp}(\alpha \mathrm{log}(\omega ))=\omega ^\alpha $$
(80)
and hence

$$\rho (\omega )\omega ^\beta $$
(81)
with $`\beta =\alpha 1`$. From an extension of this approach, one also findsComtet  that the localisation length $`\xi `$ at low frequency varies with the staggering $`V(x)`$ as $`\xi V(x)^1`$.

### VI.3 Calculation of the density of states

A central difficulty of our problem, of course, is that $`V(x)`$ does not have a simple, given distribution. Instead, it should be determined by solving Eq. (70), or calculated from the partial energy using Eq. (74), after in either case first finding the ground state configuration $`\varphi _0(x)`$, using Eq. (67) or Eq. (73). Our detailed calculations are based on Eqns. (70) and (67). Before presenting these calculations, it is useful to develop some qualitative understanding by an alternative route, using the partial energy, $`(x,\varphi )`$, and its connection to Burgers turbulence.

The evolution of $`(x,\varphi )`$ with $`x`$ is described by Eq. (72), which is similar in form to the much-studied KPZ equation.KPZ  In this correspondence, $`x`$ and $`\varphi `$ play the roles of time and space coordinates, respectively. The stochastic KPZ equation, however, which reads

$$_x+\frac{1}{2}\left(_\varphi \right)^2D_\varphi ^2=h(\varphi ,x),$$
(82)
differs from Eq. (72) in two respects. One is in the absence of the dissipation term, with coefficient $`D`$; the other is in the correlations of the force, $`h(\varphi ,x)`$. The absence of dissipation is of limited importance, because the relationship between the solutions of Eq. (82) in the small $`D`$ limit and those of Eq. (72) is well-understood: while Eq. (72) generally has solutions with many branches, corresponding to local minima in the energy of the spin chain, by taking the solution of Eq. (82) for $`D0`$ one finds the envelope of absolute ground states as a function of the boundary spin $`\varphi (x)`$. For this reason, it is rather natural to study Eq. (82) in our context.

By contrast, the nature of force correlations is more significant: while in standard form the KPZ equation has a random force $`h(\varphi ,x)`$ that is white noise in both $`x`$ and $`\varphi `$, our interest lies with correlations that are long-ranged in $`\varphi `$, and have the form

$$h(\varphi _1,x_1)h(\varphi _2,x_2)=\stackrel{~}{h}\mathrm{cos}(\varphi _1\varphi _2)\delta (x_1x_2).$$
(83)
Such correlations have been studied previously, in the context of Burgers turbulence. where the Burgers equation, Eq. (84), (equivalent to Eq. (82) with $`D0`$) together with the correlator Eq. (83), describes motion of a one-dimensional fluid.

In this analogy, we view $`\varphi (x)`$ as the coordinate of a particle as a function of time, $`x`$, and interpret Eq. (67) as the equation of motion for the particle. Imagine a fluid of such particles, without pressure and stirred randomly with force correlations derived from Eq. (83). Let $`u(\varphi ,x)=d\varphi /dx`$ be the velocity, which satisfies the Burgers equation

$$_xu+u_\varphi u=_\varphi h(\varphi ,x),$$
(84)
where comparison with Eq. (72) shows that $`u=_\varphi `$. The chiral potential is therefore given by $`V(x)=_\varphi u(\varphi _0(x),x)`$, the velocity gradient in a one-dimensional fluid flowing in a time-dependent potential $`h(\varphi ,x)`$, calculated at a point that moves with the fluid.

From literature on the Burgers equation (see Refs. Burgers, ; Pol, ; Sinai, ), or alternatively by thinking about the ground state of a spin chain as the boundary spin $`\varphi (x)`$ is varied,Fogler  one arrives at the following picture for $`(\varphi ,x)`$. As a function of $`\varphi `$ it typically consists of a few piecewise smooth sections, which meet at cusps that are local maxima of $`(\varphi ,x)`$. These cusps in $`(\varphi ,x)`$ are negative discontinuities or shocks of the Burgers velocity field, $`u(\varphi ,x)`$, at which

$$\underset{x0}{lim}\left[u(\varphi +ϵ,x)u(\varphi ,x)\right]<0.$$
(85)
They occur at the points where the ground-state spin configuration of the half-chain changes discontinuously as $`\varphi `$ varies. With increasing system length $`x`$, they undergo an evolution in which existing cusps merge and new cusps are born, at matching average rates. If the trajectory of a particle moving with the Burgers fluid should meet a shock, the particle remains trapped by the shock for all subsequent $`x`$. From these statements it is clear that (for almost all boundary conditions $`\varphi (L)`$) the ground state trajectory $`\varphi _0(x)`$ does not intersect any shocks. As a next step, from this we expect that $`_\varphi ^2(\varphi _0(x),x)>0`$, on the basis that, first, an average of $`_\varphi ^2(\varphi ,x)_\varphi u(\varphi ,x)`$ over $`\mathrm{𝑎𝑙𝑙}`$ $`\varphi `$ must be zero, since $`u(\varphi ,x)`$ is periodic in $`\varphi `$, while, second, an average restricted to $`\varphi _0`$ avoids negative discontinuities of $`u(\varphi ,x)`$. Finally, returning to excitations of the spin chain viewed as a chiral problem, we conclude that this is in the staggered regime and expect a finite localization length $`\xi _\varphi ^2(\varphi _0(x),x)>0`$; from dimensional analysis, we expect $`\xi \stackrel{~}{h}^{1/3}`$.

Two weaknesses of the argument we have given are clear: first, a more detailed treatment of averages over ground states $`\varphi _0(x)`$ would be desirable; and second, it is not certain that behavior known for chiral problems with disorder uncorrelated in $`x`$ will necessarily be present in our system, with correlations of $`V(x)`$ determined from ground-state properties. Nevertheless, results of the detailed calculations below are consistent with the foregoing conclusions.

To make further progress, we need to study statistical properties of $`V(x)`$. We find that a direct attack on this problem using the boundary conditions that are physically appropriate (in which the values of $`\varphi _0(x)`$ are specified at $`x=0`$ and $`x=L`$) is too difficult mathematically. Instead, we approach it indirectly, by relating it to a version of the problem with simpler boundary conditions, in which $`\varphi _0(x)`$ and $`_x\varphi _0(x)`$ are specified at $`x=0`$. With the latter boundary conditions, we solve jointly Eq. (67) for $`\varphi _0(x)`$ and Eq. (70) for $`V(x)`$. To discuss the connection between the two alternative sets of boundary conditions, consider Eq. (70) with an arbitrary forcing term, written as $`f(x)`$:

$$\frac{dV(x)}{dx}+V^2(x)=f(x).$$
(86)
This equation should be integrated from $`x=0`$ towards $`x=L`$, with initial condition $`V(x)\mathrm{}`$ for $`x0`$ since this is the behavior of $`_\varphi ^2(\varphi _0(x),x)`$ at small $`x`$ and $`V(x)`$ is related to the partial energy by Eq. (74). Crucially, Eq. (86) is unstable in the sense that, if $`V(x)`$ reaches sufficiently large negative values for $`V^2(x)`$ to dominate over $`f(x)`$, the solution escapes to towards $`V(x)=\mathrm{}`$, with the form $`V(x)(xx_0)^1`$ as $`x`$ approaches $`x_0`$, the position of the instability, from below. If such an instability is reached, it signals the fact that the trajectory $`\varphi _0(x)`$ no longer represents the absolute ground state of the spin chain but is instead a local maximum in the energy as a functional of configuration. If the forcing term is derived from the ground state configuration using $`f(x)=_\varphi ^2h(\varphi _0(x),x)`$, it will have correlations that ensure this instability is never reached. But if we treat Eq. (86) as a Markov process in the way we set out below, some realizations will prove unstable in this sense. Such trajectories should be discarded, and this can be arranged by supplementing Eq. (86) with absorbing boundary conditions at $`V(x)=\mathrm{}`$. The surviving trajectories must be weighted in order to sample ground states appropriately, and we return to this aspect in due course.

First we check that, specifying $`\varphi _0(x)`$ and $`_x\varphi _0(x)`$ at $`x=0`$ and integrating Eqns. (67) and (70) together, we indeed have a Markov process. This is seen most clearly by returning to the discrete equations, taking Eq. (56) in place of Eq. (70) and setting $`_\varphi ^2I1`$ to obtain

$$V_i=\frac{V_{i1}}{1+V_{i1}}+_\varphi ^2h_i(\varphi _i).$$
(87)
Similarly, in place of Eq. (73), Eq. (54) can be written as

$$\varphi _i=\varphi _{i1}+_\varphi _{i1}(\varphi _{i1}).$$
(88)
It is now clear that, since $`\varphi _i`$ is determined by $`h_j`$ for $`j<i`$, it is independent of the function $`h_i(\varphi )`$. Moreover, by construction, each $`h_j(\varphi )`$ is an independent random function. Returning to the continuum limit and noting the correlator for $`h(\varphi ,x)`$ given in Eq. (83), we see that in Eq. (86) we should take $`f(x)`$ Gaussian distributed, with zero mean and correlator

$$f(x)f(y)=\stackrel{~}{h}\delta (xy).$$
(89)
In light of the discussion given in Sec. VI.2, our next step is to find the probability of an unusually large negative fluctuation in $`V(x)`$, integrated over an interval of length $`l=x_2x_1`$. To this end, it is helpful first to define $`S(\alpha )`$ by the equation

$$\mathrm{exp}\left(\alpha _{x_1}^{x_2}𝑑xV(x)\right)=\mathrm{exp}\left(lS(\alpha )\right)$$
(90)
where the angular brackets indicate averaging over $`V`$, and we will find that $`S(\alpha )`$ is independent of $`l`$ when $`l`$ is large. We denote this probability (which appeared previously in Eq. (80)) by $`P_0(\omega )`$ below, and calculate it using

$$P_0(\omega )=\underset{l}{\mathrm{max}}_i\mathrm{}^i\mathrm{}𝑑\alpha \mathrm{exp}\left(lS(\alpha )+\alpha \mathrm{log}(\omega )\right).$$
(91)
At small $`\omega `$ we can approximate the integral by its value at the saddle-point, determined from

$$l\frac{S}{\alpha }+\mathrm{log}(\omega )=0.$$
(92)
In turn, the maximization on $`l`$ (remembering that the value of $`\alpha `$ at the saddle-point is itself a function of $`l`$) gives

$$S(\alpha )=0.$$
(93)
Therefore

$$P_0(\omega )\omega ^{\alpha _0},$$
(94)
where $`\alpha _0`$ is the solution to the Eq. (93).

Now we must compute $`S(\alpha )`$. A convenient method is to derive from the Langevin equation Eq. (86) a Fokker-Planck equation, make a similarity transformation of the latter to obtain a Schrödinger equation, and express this as a path integral. The absorbing boundary condition at $`V=\mathrm{}`$ is built in automatically in this approach. These methods are described, for example in Ref. ZinnJustin, . In this way we find

$$\mathrm{exp}\left(lS(\alpha )\right)=\frac{𝒟V\mathrm{exp}\left[_0^l𝑑x\left\{\frac{1}{2\stackrel{~}{h}}\left(\left(\frac{dV}{dx}\right)^2+V^4\right)V\left(1\alpha \right)\right\}\right]}{𝒟V\mathrm{exp}\left[_0^l𝑑x\left\{\frac{1}{2\stackrel{~}{h}}\left(\left(\frac{dV}{dx}\right)^2+V^4\right)V\right\}\right]}.$$
(95)
The path integrals in this expression are propagators in imaginary time $`x`$ for a particle moving with coordinate $`V(x)`$ in a polynomial potential, which is $`V^4V(1\alpha )`$ in the case of the numerator and $`V^4V`$ in the case of the denominator. For large $`l`$, both path integrals are dominated by ground-state contributions, justifying the $`l`$-dependence displayed in Eq. (90). In this limit, moreover, $`S(\alpha )`$ is given by the difference in ground-state energies for the two potentials. It is clear that $`S(\alpha )=S(2\alpha )`$ and that $`S(0)=0`$. Hence $`\alpha _0=2`$ and

$$P_0(\omega )\omega ^2.$$
(96)
It is incidentally also apparent that $`V>0`$, confirming our earlier argument that the chiral description of spin-chain excitations is in the staggered regime.

To complete our calculation of the density of excitations in frequency, a further step is necessary. We have so far considered spin configurations $`\varphi _0(x)`$ that are generated in the ensemble of disorder realizations using fixed values for $`\varphi _0(x)`$ and $`_x\varphi _0(x)`$ at $`x=0`$, and are local minima of the energy but not necessarily the absolute minimum. We should weight these configurations by a factor $`P_1(\omega )`$, according to the probability that they appear in an ensemble with physically appropriate boundary conditions, in which $`\varphi _0(x)`$ is fixed at $`x=0`$ and $`x=L`$. In addition, to obtain properties of excitations about the ground state, we require a further weighting factor $`P_2(\omega )`$, involving the probability that a configuration is the absolute minimum in energy. This is in principle a difficult quantity to determine,note  because it involves global features of the system, and we content ourselves with a heuristic approach which is in a similar spirit to our discussion of an anharmonic oscillator in Sec. IV.2. Specifically, we discard all configurations that have a nearby maximum of the energy. In this way we correctly exclude local minima that are separated from the ground state by the nearby maximum, but we make errors of two kinds. First, we fail to exclude local minima that have no nearby maximum and are separated by a large distance in configuration space from the true ground state: we assume (as in Sec. IV.2) that excitations about such local minima have the same statistical properties as those about the true ground state. Second, we wrongly exclude true ground states that have nearby maxima: we assume that those ground states which remain are characteristic of the whole set. We remark finally that if the factor $`P_2(\omega )`$ is omitted, we obtain an excitation density averaged over all configurations that are local energy minima.

To find the $`\omega `$-dependence of these two weighting factors, we consider a family of nearby configurations $`\varphi (s,x)=\varphi _0(x)+\eta (s,x)`$ for $`x_1xx_2`$, parameterized by $`s`$. The corresponding family of chiral potentials is $`V(s,x)=V(x)+W(s,x)`$, and we restrict attention to small $`\eta (s,x)`$ and $`W(s,x)`$. The weight $`P_1(\omega )`$ appears because, in disorder realizations which generate negative fluctuations of $`V(x)`$, trajectories of $`\varphi (s,x)`$ as a function of $`x`$ are compressed by an amount which increases with the size of the integrated potential fluctuation. It is therefore determined by comparing $`\eta (s,x_2)`$ with $`\eta (s,x_1)`$. The weight $`P_2(\omega )`$ is determined by finding the probability for escape of $`W(s,x)`$ to negative infinity, signaling the occurrence of an energy maximum.

We find the evolution of $`\eta (s,x)`$ and $`W(s,x)`$ in terms of $`V(x)`$ by linearizing Eqns. (67) and (70) respectively, obtaining

$$\frac{d\eta }{dx}=V(x)\eta (s,x)$$
(97)
and

$$\frac{dW}{dx}+2V(x)W(s,x)=_\varphi ^3h(\varphi ,x)\eta (s,x).$$
(98)
These equations have the solutions

$$\eta (s,x)=\mathrm{exp}\left[_{x_1}^x𝑑yV(y)\right]\eta (s,x_1)$$
(99)
and

$$W(s,x)=\mathrm{exp}\left(2_{x_1}^x𝑑yV(y)\right)\left(W(s,x_1)+\eta (s,x_1)_{x_1}^x𝑑y\left[\mathrm{exp}\left(3_{x_1}^z𝑑zV(z)\right)_\varphi ^3h(\varphi _0(y),y)\right]\right).$$
(100)
Let us choose the parametrization $`s`$ in such a way that at $`s=0`$, $`\eta (s,x_1)=0`$ and $`V(s,x_1)`$ has a minimum. Then for small $`s`$, $`\eta (s,x_1)s`$ and $`W(s,x_1)s^2`$. We find $`P_1(\omega )`$ as follows. For the integrated value of $`V(s,x)`$ along a trajectory from the family to have a value similar to that at $`s=0`$, the integral of $`W(s,x_2)`$ must not be too large. We hence select those values of $`s`$ for which $`_{x_1}^{x_2}𝑑xW(s,x)1`$. Using Eq. (100) and that fact that $`_{x_1}^{x_2}𝑑xV(0,x)=\mathrm{log}(\omega )`$, this requires $`|s|\omega `$, which in turn implies that $`|\eta (s,x_1)|\omega `$ and $`|\eta (s,x_2)|\omega ^2`$. Since the physically relevant boundary conditions for the spin chain fix the value of $`\varphi (x)`$ at $`x=L`$ (that is, at larger $`x`$), we weight configurations uniformly in $`\eta (s,x_2)`$ and conclude that $`P_1(\omega )\omega ^2`$. Turning to $`P_2(\omega )`$, we require that $`W(s,x_2)`$ calculated using our linearization should not be large and negative for any $`s`$, since otherwise a full treatment, including non-linearities, would with high probability result in escape. Minimizing $`W(s,x_2)`$ with respect to $`s`$, using Eq. (100), we require

$$\underset{s}{\mathrm{min}}W(s,x_2)\omega ^2I^21,$$
(101)
where

$$I=_{x_1}^{x_2}𝑑x\left[\mathrm{exp}\left(3_{x_1}^x𝑑yV(y)\right)_\varphi ^3h(\varphi _0(x),x)\right].$$
(102)
For small $`\omega `$ and large $`|x_2x_1|`$, $`I`$ is a random variable whose distribution is independent of $`\omega `$, and hence the probability that $`\mathrm{min}_sW(s,x_2)1`$ is $`P_2(\omega )\omega `$.

Combining these results with the form derived for $`P_0(\omega )`$, we find an integrated density of states for excitations, averaging over all local minima of Eq. (46), which varies as

$$N(\omega )P_0(\omega )P_1(\omega )=\omega ^4,$$
(103)
while for the global minimum we find

$$N(\omega )P_0(\omega )P_1(\omega )P_2(\omega )=\omega ^5.$$
(104)
Hence $`\beta =3`$ or $`\beta =4`$, in agreement with Eqs. (43) and (45). Strikingly, the behavior derived here for the random field $`XY`$ chain matches that expected from the simple discussion of an anharmonic oscillator, given in Sec. IV.2.

## VII Other systems without Goldstone modes

In this section we discuss the extent to which the approach we have set out for the random field $`XY`$ spin chain can be extended to treat excitations in other models for disordered systems without a broken continuous symmetry. The $`XY`$ spin chain has two obvious and important simplifying features: disorder couples only to one of the dynamical variables ($`\varphi `$ but not $`\mathrm{\Pi }`$), so that frequencies are simply related to stiffnesses; and the system is one-dimensional. The calculations we have described involved several steps: mapping to a chiral formulation; analytic determination of the excitation density in frequency via a study of statistical properties of disorder within this chiral formulation; and, for an efficient numerical treatment, the use of a partial energy. As we broaden the range of models under consideration, fewer of these steps remain possible. We consider, first, a generic one-dimensional continuum system, with a pair of conjugate dynamical variables at each point and disorder that couples to both, and second, a specific example of such a problem, the random field Heisenberg spin chain. In these cases we find a chiral description and introduce a partial energy, but are not able to determine statistical properties of the disorder appearing in this description. Instead, we derive for such problems the general relation between densities of stiffnesses and of excitation frequencies, suggested from an analysis of random matrix theory in Eq. (22) above. Third, for the random field $`XY`$ model in two dimensions we show how to introduce a chiral description, but leave applications of this for future work. For all these problems, we expect $`\rho (\omega )\omega ^4`$ at small $`\omega `$, on the general grounds discussed in Sec. III.

### VII.1 Generic one-dimensional problem

Consider a one-dimensional system with conjugate dynamical variables $`\varphi _1(x)`$ and $`\varphi _2(x)`$, and the Hamiltonian

$$H=𝑑x\left[\frac{1}{2}\left\{\left(\frac{d\varphi _1}{dx}\right)^2+\left(\frac{d\varphi _2}{dx}\right)^2\right\}+h(\varphi _1,\varphi _2,x)\right].$$
(105)
We wish to study harmonic excitations about the ground state: for $`i=1,2`$, let $`\varphi _i^0`$ be the configuration that minimises the energy, Eq. (105), and write $`\varphi _i=\varphi _i^0+\psi _i`$. Expanding to quadratic order in $`\psi _i`$,

$$=\left(\begin{array}{cc}_x^2+_{\varphi _1}^2h& _{\varphi _1}_{\varphi _2}h\\ _{\varphi _1}_{\varphi _2}h& _x^2+_{\varphi _2}^2h\end{array}\right).$$
(106)
To write this as $`=Q^TQ`$ and find $`Q`$, we essentially repeat the sequence of arguments which led from Eq. (66) to Eq. (75). We introduce a partial energy $`(\varphi _1,\varphi _2,x)`$ which satisfies

$$_x+\frac{1}{2}\left\{\left(\frac{}{\varphi _1}\right)^2+\left(\frac{}{\varphi _2}\right)^2\right\}=h(\varphi _1,\varphi _2,x)$$
(107)
and may be used to find the ground state configuration via the analog of Eq. (73),

$$\frac{d\varphi _i^0}{dx}=_{\varphi _i}E(\varphi _1^0,\varphi _2^0,x).$$
(108)
Then we define a $`2\times 2`$ matrix version of the chiral potential, via

$$V_{ij}=_{\varphi _i}_{\varphi _j}(\varphi _1^0,\varphi _2^0,x).$$
(109)
It satisfies

$$\frac{d}{dx}V_{ij}+\underset{k=1}{\overset{2}{}}V_{ik}V_{kj}=_{\varphi _i}_{\varphi _j}h(\varphi _1,\varphi _2,x)$$
(110)
and therefore $`Q`$ can be chosen to have elements

$$Q_{ij}=\frac{d}{dx}\delta _{ij}+V_{ij}.$$
(111)
In this way, we have written Eq. (106) as a type of chiral problem, which we term two-channel because $`V`$ is a $`2\times 2`$ matrix.

Further analysis can be divided into two stages. One stage, given $`V`$, is to find $`d(\kappa )`$ and $`\rho (\omega )`$, as was done for the single-channel problem in Sec. VI.2. The other stage is to determine the distribution for $`V`$, as was done for the single-channel problem in Sec. VI.3. General one-dimensional multichannel chiral problems of the type that arise in our calculation of stiffnesses have been studied in Refs. Brouwer,  and Brouwer2, , with $`V`$ chosen Gaussian distributed and uncorrelated in $`x`$. However, as far as we are aware, there has been no previous work on multichannel chiral problems of the type that give frequencies. For the two-channel problem, while we have not been able to make progress in obtaining the distribution of $`V`$, we have been able to find a general relation between calculations of stiffnesses and frequencies. This connects $`d(\kappa )`$ and $`\rho (\omega )`$ in the way given by Eq. (22). We present these arguments next.

### VII.2 Two-channel Bosonic Chiral Problems

In this subsection we compare the calculation of stiffnesses $`\kappa \lambda ^2`$, determined from the eigenvalue problem

$$\left(\begin{array}{cc}0& \frac{d}{dx}+V\\ \frac{d}{dx}+V^T& 0\end{array}\right)\left(\begin{array}{c}\psi _1\\ \psi _2\end{array}\right)=\lambda \left(\begin{array}{c}\psi _1\\ \psi _2\end{array}\right),$$
(112)
with the calculation of frequencies $`\omega ϵ^2`$, determined from

$$\left(\begin{array}{cc}0& \frac{d}{dx}+V\\ \frac{d}{dx}+V^T& 0\end{array}\right)\left(\begin{array}{c}\psi _1\\ \psi _2\end{array}\right)=ϵ\left(\begin{array}{c}\psi _1\\ \sigma _2\psi _2\end{array}\right).$$
(113)
We do this following the technique set out for the single-channel problem after Eq. (76), and described for multichannel problems in Ref. Brouwer2, . For both stiffnesses and frequencies, we write $`\psi _1=a\psi _2`$, where $`a`$ is an $`x`$-dependent $`2\times 2`$ matrix. The eigenvalue density, $`d(\kappa )`$ or $`\rho (\omega )`$, is then determined from the evolution of $`a`$ with $`x`$, via a node-counting theorem. The evolution equation for $`a`$ is, in the case of stiffnesses, from Eq. (112),

$$\frac{da_s}{dx}=\lambda \left(𝟙+𝕒_𝕤^\mathrm{𝟚}\right)a_sVV^Ta_s,$$
(114)
and in the case of frequencies, from Eq. (113),

$$\frac{da_f}{dx}=ϵ\left(\sigma _2+a_f^2\right)a_fVV^Ta_f.$$
(115)
These equations are the equivalent for the two-channel problem of Eq. (77) for the single-channel problem. Following Ref. Brouwer2, , the density of states is given by the rate at which the eigenvalues of $`a_s`$ or $`a_f`$ move with increasing $`x`$ along the real axis. Our aim is to compare this rate in the two cases, taking $`\kappa `$ and $`\omega `$ small, in a way that does not require detailed knowledge of the distribution of $`V`$. We assume only that (as for single-channel problem) $`V`$ fluctuates about a non-zero mean, and that both the fluctuations and the mean of $`V`$ have comparable importance in the evolution of $`a_s`$ and $`a_f`$ with $`x`$.

In outline, this evolution is as follows. If fluctuations in $`V`$ are omitted, $`a_s`$ and $`a_f`$ have stable fixed points at which both their eigenvalues are small in $`\kappa `$ and $`\omega `$ respectively. In both cases, the stable fixed point has a basin of attraction with a boundary that is reached when an eigenvalue of $`a_s`$ or $`a_f`$ is large and positive ($`𝒪(\lambda ^1)`$ or $`𝒪(ϵ^1)`$ respectively). Restoring fluctuations in $`V`$, with increasing $`x`$ we find (partly on the basis of simulations, not presented here) that $`a_s`$ or $`a_f`$ fluctuates in the vicinity of its fixed point for intervals that are long if $`\kappa `$ or $`\omega `$ are small. Each such interval ends when fluctuations take the matrix to the boundary of the basin of attraction. One eigenvalue of $`a_s`$ or $`a_f`$ then runs off to positive infinity, reappears from negative infinity and returns to the vicinity of the fixed point. This is entirely analogous to the evolution in the single-channel problem of $`\mathrm{tan}(\theta )`$ as $`\theta `$ increases from $`\theta n\pi `$ to $`\theta (n+1)\pi `$, which is described in Sec. VI.2.

To develop the picture further, we consider separately the regions close to the fixed points, which are different for $`a_s`$ and $`a_f`$, and the region far from the fixed points, which is essentially the same in both cases. It is useful to introduce an explicit coordinate system. We write

$$V=\frac{1}{2}(𝟙+\xi _\mathrm{𝟘}𝟙+\xi _\mathrm{𝟙}\sigma _\mathrm{𝟙}+\xi _\mathrm{𝟛}\sigma _\mathrm{𝟛}),$$
(116)
where $`\xi _0`$, $`\xi _1`$ and $`\xi _3`$ are random with mean zero, and $`V`$ is taken proportional to the unit matrix, with a proportionality constant that can be changed by a rescaling of the length coordinate $`x`$. We also set

$$a_{s,f}=s𝟙+𝕫_\mathrm{𝟙}\sigma _\mathrm{𝟙}+𝕫_\mathrm{𝟚}\sigma _\mathrm{𝟚}+𝕫_\mathrm{𝟛}\sigma _\mathrm{𝟛}.$$
(117)
For calculation of stiffnesses using this coordinate system, one has $`z_2=0`$ and Eq. (114) becomes

$`{\displaystyle \frac{dz_1}{dx}}=2\lambda sz_1z_1+\xi _0z_1+\xi _1s,`$ (118)
$`{\displaystyle \frac{dz_3}{dx}}=2\lambda sz_3z_3+\xi _0z_3+\xi _3s,`$
$`{\displaystyle \frac{ds}{dx}}=\lambda \left(s^2+z_1^2+z_3^2\right)+\lambda s+\xi _0s+\xi _1z_1+\xi _3z_3.`$
If fluctuations in $`V`$ are omitted (by setting $`\xi _i=0`$ for $`i=0,1`$ and $`3`$), these equations have a stable fixed point at $`s\lambda `$, $`z_1=z_3=0`$. Including fluctuations in $`V`$, the typical magnitudes of $`s`$, $`z_1`$ and $`z_3`$ are $`𝒪(\lambda )`$.

By contrast, for calculation of frequencies Eq. (115) gives

$`{\displaystyle \frac{dz_1}{dx}}=2ϵsz_1z_1+\xi _0z_1+\xi _1s,`$ (119)
$`{\displaystyle \frac{dz_2}{dx}}=2ϵsz_2+ϵz_2+\xi _0z_2,`$
$`{\displaystyle \frac{dz_3}{dx}}=2ϵsz_3z_3+\xi _0z_3+\xi _3s,`$
$`{\displaystyle \frac{ds}{dx}}=ϵ\left(s^2+z_1^2+z_2^2+z_3^2\right)s+\xi _0s+\xi _1z_1+\xi _3z_3.`$
In this case, if fluctuations in $`V`$ are omitted there is a fixed point at $`sϵ^3`$, $`z_2ϵ`$, $`z_1=z_3=0`$. Including fluctuations, the typical magnitudes of $`s`$, $`z_1`$ and $`z_3`$ are $`𝒪(ϵ^3)`$, while $`z_2`$ remains $`𝒪(ϵ)`$.

Now consider escape of $`a_s`$ or $`a_f`$ from the vicinity of the relevant fixed point, which requires $`s\lambda `$ in the case of stiffnesses and $`sϵ^3`$ in the case of frequencies. In these regimes we argue that the evolution equations in the two instances are essentially equivalent. More specifically, we make two approximations. First, we take $`z_2=0`$ in both cases, even though this is exact only in the first case. We do so because in the second case it is clear from Eq. (119) that $`z_2ϵ`$ provided $`s<(2ϵ)^1`$, and hence that non-zero $`z_2`$ has no important effect on the evolution of $`z_1`$, $`z_3`$ and $`s`$. Second, for $`a_s`$ or $`a_f`$ far from its fixed point, we omit the terms independent of $`a_s`$ or $`a_f`$ (and small in $`\kappa `$ or $`\omega `$) from the right-hand sides of Eqns. (114) and (115). With these approximations, and making the rescalings $`a=\lambda a_s`$ and $`a=ϵa_f`$, the two equations both become

$$\frac{da}{dt}=a^2aVV^Ta.$$
(120)
This stochastic process results in a stationary probability distribution for $`a`$ if absorbing boundary conditions are imposed at infinity and probability flux is injected near $`a=0`$. In the stationary state there is an eigenvalue flux along the positive real axis, with a rate that is determined by the probability density in a region near $`a=0`$. The size of this region, and therefore the probability density within it, follows in each case from our discussion of the fixed points and their neighborhoods. Characterizing $`a`$ by the value of its trace, and combining the fixed-point coordinates with the rescalings used to arrive at Eq. (120), the value of $`\mathrm{Tr}a`$ in this region is $`𝒪(\lambda ^2)`$ in the first instance and $`𝒪(ϵ^4)`$ in the second. From this, we conclude that the integrated densities of stiffnesses and frequencies are equal for $`\lambda ^2ϵ^4`$. We hence recover from this discussion of one-dimensional systems the result given in Eq. (22), which was reached in Sec. III by a very different route.

For completeness, we remark that the derivation presented here crucially depends on $`V(x)0`$, and assumes that $`\xi _0`$, $`\xi _1`$, and $`\xi _3`$ become comparable to 1 only during their rare large fluctuations. And indeed, for $`V(x)=0`$ it has been conjectured and checked numericallyGAunp  that in that case $`\rho (\omega )\omega ^{1/3}`$ even though it is well knownBrouwer2  that $`d(\kappa )\mathrm{log}(\kappa )`$.

### VII.3 The random field Heisenberg spin chain

The random field Heisenberg spin chain is perhaps the most obvious example of the generic one-dimensional problems discussed above. The conjugate dynamical variables are, of course, the two components for displacement of a spin from its orientation in the ground state, and phase space is made up of spheres for each point on the chain. The curvature of phase space is responsible for some changes in formuls derived in Sec. VII.1, which we now set out.

The continuum version of the random field Heisenberg spin chain, parameterizing spin orientations using the angles $`\theta (x)`$ and $`\varphi (x)`$, is

$$H=𝑑x\left[\frac{1}{2}\left\{\left(\frac{d\theta }{dx}\right)^2+\mathrm{sin}^2\left(\theta \right)\left(\frac{d\varphi }{dx}\right)^2\right\}+h(\varphi ,\theta ,x)\right].$$

The evolution equation for the partial energy is

$$\frac{}{x}+\frac{1}{2}\left\{\left(\frac{}{\theta }\right)^2+\frac{1}{\mathrm{sin}^2(\theta )}\left(\frac{}{\varphi }\right)^2\right\}=h(\varphi ,\theta ,x)$$

and, given $``$, the ground state configuration can be calculated using

$$\frac{d\varphi _0}{dx}=\frac{1}{\mathrm{sin}^2(\theta )}_\varphi ,\frac{d\theta _0}{dx}=_\theta .$$

Introducing coordinates for deviations from the ground-state configuration, $`\theta =\theta _0+\psi _\theta `$ and $`\varphi =\varphi _0+\psi _\varphi `$, we take as canonical variables $`\psi _\theta `$ and $`\psi _\varphi \mathrm{sin}\left(\theta _0\right)`$. The quadratic Hamiltonian can be written $`=Q^TQ`$, with $`Q`$ given by Eq. (111) in terms of the $`2\times 2`$ matrix $`V`$, which can be computed from $``$ using

$$V=\left[\begin{array}{cc}\mathrm{cot}(\theta _0)+\frac{1}{\mathrm{sin}^2(\theta _0)}_\varphi ^2& \frac{1}{\mathrm{sin}^2(\theta _0)}_\theta _\varphi 2\frac{\mathrm{cos}(\theta _0)}{\mathrm{sin}^3(\theta _0)}_\varphi \\ \frac{1}{\mathrm{sin}(\theta _0)}_\varphi _\theta & _\theta ^2\end{array}\right].$$

This formalism would provide the starting point for a treatment of the one-dimensional random field Heisenberg model similar to that presented for the $`XY`$ model in Sec. V.

### VII.4 Two-dimensional random field $`XY`$ model

The two-dimensional random field $`XY`$ model provides an example for which our formulation of excitation problems using chiral Hamiltonians can be carried through in more than one dimension. The two-dimensional version of Eq. (66) is

$$K=\left\{\frac{1}{2}\left[\mathrm{\Pi }^2+\left(_x\varphi \right)^2+\left(_y\varphi \right)^2\right]+h(\varphi ,x,y)\right\}𝑑x𝑑y$$

and the ground state $`\varphi _0`$ satisfies

$$\left[\frac{^2}{x^2}+\frac{^2}{y^2}\right]\varphi _0+_\varphi h(\varphi _0,x,y)=0.$$
(121)
The amplitudes of normal mode excitations obey

$$\left[\frac{^2}{x^2}\frac{^2}{y^2}+_\varphi ^2h(\varphi ,x,y)\right]\psi =E\psi .$$
(122)
We introduce a chiral potential $`V(x,y)`$, defined as the solution to

$$_x^2V+(_xV)^2+_y^2V+(_yV)^2=_\varphi ^2h(\varphi _0,x,y),$$
(123)
and the matrix $`Q`$, given by

$$Q=\left(\begin{array}{cc}_x+_xV& _y_yV\\ _y+_yV& _x+_xV\end{array}\right).$$
(124)
Then the matrix $`Q^TQ`$ is diagonal, with the operator of Eq. (121) as diagonal elements. The construction of a chiral matrix $`\stackrel{~}{}`$ is therefore complete for this problem, too: $`\stackrel{~}{}`$ in fact has the form of a random Dirac Hamiltonian, studied in Ref. LGL, . While we have not been able to find a quantity for the two-dimensional system that plays the role of the partial energy in one dimension, we have formulated a generalization of the Burgers equation, which is sufficient to show that a solution to Eq. (123) exists. The practical determination of the chiral potential, however, remains an open problem in this case.

## VIII Relation to experiment

The range of physical systems to which the ideas we have set out may apply is potentially very wide. The essential requirements are: quenched disorder; a ground state that can be considered classically; and excitations that can be treated as weakly interacting normal modes. These requirements may be met for vibrational modes, either in glasses, or in randomly pinned phases with broken translational symmetry such as charge density wave states. They may also be satisfied for magnetic excitations in disordered ferromagnets or antiferromagnets, or in spin glasses, either with or without significant magnetic anisotropy. The most specific feature for which one would like experimental evidence is probably the frequency dependence $`\rho (\omega )\omega ^4`$ of the density for excitations that are not Goldstone modes. Alternatively, for excitations that are Goldstone modes, interest focusses on the frequency dependence of the mean free path at small $`\omega `$, or, in the case of Heisenberg antiferromagnets and spin glasses which are quasi one-dimensional, the low-frequency form of the density of states.

Before summarising the experimental situation, it is useful to outline ways in which excitations in a disordered system may fall outside the framework we have used. Large quantum fluctuations are the obvious route to different physics, and may be important in at least two ways. First, it can happen that a ground state is very far from being a classical configuration dressed with small zero-point fluctuations: random singlet phases SH  in disordered quantum spin chains constitute a well-studied example. Second, it may be that even low-lying quantized excitations are unlike weakly interacting bosonic modes. To have a concrete example, consider the single mode problem for an anharmonic oscillator with potential $`U(q)`$, as reviewed in Sec. IV.2. Typically, this potential energy will have local minima separated by barriers from the absolute minimum, and the importance of quantum motion through or over the barrier will depend on the size of the mass $`m`$ appearing in Eq. (38). Thus, for a given $`U(q)`$, in the semi-classical limit the low-lying quantum states will be close to harmonic oscillator levels with the classical normal mode frequency, while if quantum fluctuations are large, tunneling through barriers may hybridise levels associated with different classical minima, generating two-level systems, or the low-lying levels may be determined by the form of $`U(q)`$ in regions too far from its absolute minimum for the classical normal mode frequency to be relevant. In the ensemble, this crossover takes place as a function of frequency, as discussed in Ref. IK, . As a result, one expects harmonic excitations with a density $`\rho (\omega )\omega ^4`$ at higher frequencies, and two-level systems with a constant density at lower frequencies.

Turning to experiments, somewhat surprisingly, the best evidence that we are aware of for an excitation density with an $`\omega ^4`$ dependence is from studies of vibrational excitations in glasses. Here, inelastic neutron scattering and other measurements BuchenauPRL  indicate an excess density of harmonic modes, compared to what is expected from Debye theory and the measured speed of sound. Analysis of the temperature dependence of the heat capacity Buchenau  separates two contributions that are additional to the Debye value. One is approximately linear in temperature $`T`$, dominates at low $`T`$, and is attributed to two-level systems. The other varies as $`T^5`$, dominates at higher $`T`$, and is attributed to harmonic modes with the stated frequency dependence for their density. From a theoretical viewpoint, such excitations differ in an important way from the one we have discussed, because they coexist with propagating phonon modes. The consequences of coupling between the two sets of modes remain to be investigated, although coupling between localised anharmonic vibrational modes and phonons has been studied in Ref. Yu, .

A context in which localised vibrational modes are expected without any coexisting propagating phonons is provided by pinned charge density waves, represented in one dimension by the model studied in Sec. VI. Indeed, it was in this framework that Aleiner and Ruzin AR  and Fogler Fogler  argued for a density of harmonic proportional to $`\omega ^4`$, with a crossover at low frequency Fogler  to a constant density of two-level systems. As reviewed by Fogler Fogler , existing measurements of frequency-dependent conductivity Wu  do not show the response expected from such harmonic modes, possibly because the low-temperature limit is not accessed.

Studies of spin waves in disordered magnetic systems present opportunities to examine both Goldstone and non-Goldstone modes. In particular, inelastic neutron scattering measurements of spin dynamics in a dilute near-Heisenberg antiferromagnet Uemura  show the expected broadening in wavevector of excitations with increasing frequency. In contrast, magnetic neutron scattering measurements on amorphous magnetic alloys in which there is local magnetic anisotropy Cowley  find modes which are broad in wavevector at all frequencies. In this case the density of excitations is approximately constant in frequency over the measured range. From our results, we expect a decrease in this density at low frequency, and a more detailed examination of low-frequency behavior would be of considerable interest.

## IX Concluding remarks

Since we have investigated a variety of different directions in this paper, it is perhaps useful to close with a short summary of our main points.

Considering the general quadratic Hamiltonian for bosonic excitations, Eq. (7), we have noted a formal similarity between it and fermionic random matrix Hamiltonians from the additional symmetry classes (e.g. Eq. (5)). We have argued that this has limited direct consequences, because stability of the ground state requires correlations between matrix elements of the bosonic Hamiltonian which invalidate a random matrix approach. Instead we have shown that there is a useful mapping to an auxiliary problem with the structure of the chiral symmetry class, which we have set out explicitly for a range of models.

Examples of bosonic excitations in disordered media separate into those that are Goldstone modes, and those that are not. For the former, we have reviewed established results, which demonstrate that low frequency excitations decouple from disorder except for some systems in and below a marginal dimension $`d_c=2`$. For excitations that are not Goldstone modes, we have underlined the way in which disorder itself generates low frequency excitations, and the universal form expected for their density, $`\rho (\omega )\omega ^4`$. Taking as an example the one-dimensional $`XY`$ model in a random field, we have used our techniques in detailed analytical and numerical calculations, obtaining results that illustrate this behavior of $`\rho (\omega )`$.

A further application of the techniques we have developed here is to Mattis models. These model spin glasses lack frustration, and their statistical mechanical properties are equivalent under a gauge transformation to those of ferromagnets. They have a ground state spin configurations which are known for each disorder configuration, but excitations are nevertheless affected in a non-trivial way by disorder.Sher  Because the ground state is known explicitly, the magnons in Mattis glasses are much easier to study than in a real spin glass. Being Goldstone modes, these magnons fall into the same category as those in weakly disordered antiferromagnets, discussed in this paper. In particular, their critical dimension is $`d_c=2`$. In future work AG  by one of the present authors and A. Altland, the localization and transport properties of magnons in two and three dimensional Mattis glasses will be explored.

## Acknowledgments

We thank A. Altland and M. Zirnbauer for valuable discussions. The work was supported by EPSRC under Grant No. GR/J78327, and in part by the National Science Foundation under Grant No. PHY99-07949.