text:
# Particle Acceleration through Multiple Conversions from Charged into Neutral State and Back

## I Introduction

The cosmic rays (CRs) have been studied for a long time (see, e.g., CRrew  for a review). There are many models of their (likely non-uniform) sources, but explaining the origin of the highest-energy particles ($`10^{20}`$ eV) observed in CRs is still a challenging problem prd . In general, two scenarios have been proposed for the particle acceleration in astrophysical environments. One is the acceleration by electric field in geometries, where this field is not perpendicular to magnetic field-lines, for example, in the vicinity of magnetized rotating neutron stars NSacc . Very hard particle spectra may emerge in this way, but, because of the curvature losses inherent to such geometries, the upper limit to proton energy appears to be below $`10^{20}`$ eV.

Another class of scenarios assumes gradual, “stochastic” (or diffusive) acceleration of charged particles through multiple reflections from inhomogeneities of magnetic field in the environments where large velocity gradients are present (see fermi  for a review). According to the generally accepted view (see, e.g., fermi ; flow1 ; flow2 ) this mechanism works equally well for shocks and shear flows, and the emerging spectra of accelerated particles are such, that only a small fraction of the total energy budget is contained in the most energetic particles.

Furthermore, there is a serious obstacle in achieving high energies of the accelerated particles – their diffusive escape. In the environments, where bulk velocities are sub-relativistic ($`vc`$), it takes many reflections to increase the energy of accelerated particle twofold. In this case, the mechanism can work only in the largest objects in the Universe, like galaxy clusters or radio lobes and knots in active galactic nuclei (AGNs), and only under optimistic assumption about the diffusion coefficient for particles in magnetic field. For example, Bohm diffusion leads to particle escape on the timescale $`R_0^2/r_gc`$, where $`R_0`$ is the accelerator’s size and $`r_g`$ the particle’s gyroradius, while the acceleration occurs on the timescale $`(c/v)^2r_g/c`$. The maximum attainable energy is defined by equating these timescales: $`r_g=(v/c)R_0R_0`$. It is of the order of $`10^{19}`$ eV in the most favorable cases.

Ultrarelativistic shocks and shear flows are more promising, since the parameter $`v/c1`$. Moreover, the energy gain per cycle, i.e., at each reflection from the shock or from the shear flow boundary, can approach the factor of $`\mathrm{\Gamma }^2`$, where $`\mathrm{\Gamma }`$ is the bulk Lorentz factor (measured in the upstream fluid frame in the case of shock). Such a scenario was suggested, e.g., as the dominant acceleration mechanism in gamma-ray bursts (GRBs) LevEich ; Viet ; Waxm . However, in the standard acceleration theory, the energy gain of $`\mathrm{\Gamma }^2`$ occurs only at the first cycle, while the subsequent ones result in energy gain of $`2`$ each BO ; kirk .

The reason is that the shock catches up to the reflected particle (or the particle crosses the boundary of the shear flow) when its trajectory makes an angle of $`2/\mathrm{\Gamma }`$ to the shock normal (or bulk velocity vector), while keeping the $`\mathrm{\Gamma }^2`$ energy gain would require isotropisation.

The smallness of this angle is apparent for the case of shear flow (see Fig. 1b). For shocks it is a consequence of relativistic motion (Fig. 1a). Indeed, the component of particle velocity along the shock normal becomes smaller than the velocity of the shock itself as soon as the angle between the particle’s momentum and the shock normal grows larger than $`1/\mathrm{\Gamma }`$. There is no apparent way to isotropize reflected particles unless a special structure in the magnetic field ahead of the shock is introduced, e.g., a counter-propagating shock.

However, a means to circumvent these limitations exists Rome  – it can be done by switching particle’s charge on and off at right times. Paradoxically, interactions with photons, which have been always treated as dissipative processs leading to degradation of particle energy, in fact play positive role: they allow (through the charge-changing particle conversion) to keep the $`\mathrm{\Gamma }^2`$ energy gain up to the largest particle energies. There are also other types of conversion reactions, which we briefly discuss in the following section. One of the suggested below charge-changing schemes (electron-photon-electron reaction chain in GRB shocks) is independently considered in Stern .

In order to outline the general picture, we intentionally skip some details, which are not essential for the proposed acceleration mechanism, but may change its quantitative characteristics. In particular, we assume that the magnetic field is either chaotic (turbulent) or uniform with field lines perpendicular to both the momentum of accelerated particle and the velocity of the flow, and treat shocks and shear flows as one-dimensional discontinuities.

We use the convention $`F_\nu \nu ^q`$ to define the spectral index $`q`$, where $`F_\nu `$ is the energy flux per unit frequency interval.

## II The converter acceleration mechanism

Two basic types of photon-induced conversion involve nucleons or electrons/positrons. Both cycles consist of two reactions:

$$p+\gamma n+\pi ^+\text{and}n+\gamma p+\pi ^{}$$
(1)
for proton cycle,

$$e^\pm +\gamma e^\pm +\gamma ^{}\text{and}\gamma ^{}+\gamma e^++e^{}$$
(2)
for electron cycle. Here $`p`$, $`n`$, $`\pi ^\pm `$, and $`e^\pm `$ denote proton, neutron, charged pions (positive and negative), positron and electron, respectively; $`\gamma `$ is a low-energy background photon and $`\gamma ^{}`$ the high-energy comptonized photon.

The second of reactions from the electron cycle (2) has a kinematic threshold $`\mathrm{\Delta }_e=2m_ec^2`$ in the center-of-momentum frame, where $`m_e`$ is the electron mass. Effectively, the first reaction also has the threshold $`\mathrm{\Delta }_e`$, since at lower energies of incident photons the efficiency of energy transfer to the comptonized photon becomes much less than unity. The reactions proceed differently depending on the background photon spectrum. Soft spectrum blocks the electron cycle as the fraction of energy transferred to a comptonized photon is too small and there are few target photons sufficiently energetic for the second reaction of the electron cycle. An example of soft spectrum could be a power-law with spectral index $`q<1`$ or a narrow-band spectrum, like black-body or line emission, with typical photon energy $`\overline{\epsilon }m_e^2c^4/\epsilon _e`$, where $`\epsilon _e`$ is the electron energy. For intermediate spectra (e.g., power-laws with indices $`1<q<1`$ or narrow-band spectrum with $`\overline{\epsilon }m_e^2c^4/\epsilon _e`$), the comptonized photon takes about 1/2 of the electron (positron) energy, and in the consequent pair production event this energy is divided nearly in equal parts between the daughter electron and positron. The cross-section in both processes is $`\sigma _{e,\gamma }10^{25}`$ cm2. For hard spectra ($`q>1`$ or $`\overline{\epsilon }m_e^2c^4/\epsilon _e`$), these reactions proceed in the deep Klein-Nishina regime, i.e., the comptonization and pair-production cross-sections decrease inversely proportional to the square of the center-of-momentum energy (their ratio is 1:2) and almost all the energy of interacting particles is transferred to one of the daughter particles. In effect, the energy losses for the combined electron/photon particle become very gradual. This case is the closest to the pure conversion (the probability of charge change $`p_c=1`$) without accompanying energy losses provided the synchrotron emission is negligible.

The reactions from the proton cycle (1) (see, e.g., RM ) have kinematic threshold of $`(m_\pi c^2+m_\pi ^2c^2/2m_N)150`$ MeV in the nucleon rest frame ($`m_\pi 140`$ MeV$`/c^2`$ is the charged pion mass and $`m_N940`$ MeV$`/c^2`$ the nucleon mass). Side by side with the reactions (1) proceed other photopionic reactions with formation of neutral pion, which preserve nucleon’s charge. They have roughly the same cross-section and should be considered as a background energy losses. The total photopionic cross-section rapidly increases with energy of incident photon and reaches maximum of $`\sigma _\pi 6\times 10^{28}`$ cm2 at $`\mathrm{\Delta }_p340`$ MeV, which corresponds to formation of $`\mathrm{\Delta }`$-resonance and should be considered as effective threshold. Well above the resonance energy, the cross-section decreases and levels off at $`10^{28}`$ cm2. The probability of charge change in a photopionic reaction is $`p_c1/3`$ at the resonance and $`p_c1/2`$ at the plateau. The inelasticity is $`0.2`$ and about 0.5, respectively.

A competing photon-induced reaction is the process of creation of an electron-positron pair by a photon interacting with the electric field of a proton $`p+\gamma p+e^{}+e^+`$, which has the cross-section $`5\times 10^{27}`$ cm2 and inelasticity $`10^3`$. With the decrease of spectral index $`q`$ of target photon field, this process becomes an increasingly important energy loss channel, and at $`q1.55`$ the difference with the photopionic processes in the inelasticity and cross-section is exactly balanced by larger number of target photons (thanks to lower threshold). Anyway, the $`p+\gamma p+e^{}+e^+`$ process can be neglected for spectra with $`q1.5`$.

In dense environments and at relatively low nucleon energies, i.e., in the case where there are few target photons, the proton cycle proceeds through inelastic nucleon-nucleon collisions, for example,

$$p+pn+p+\pi ^+\text{and}n+pp+p+\pi ^{}.$$
(3)
The kinematic threshold for these reactions is $`m_\pi c^2`$ in the center-of-momentum frame, and the cross-section at energies $`m_\pi c^2`$ is $`3\times 10^{26}`$ cm2. The acceleration via inelastic nucleon-nucleon collisions could be important in GRB internal shocks, where the required column density of $`10`$ g/cm2 is achieved insh .

In both proton and electron cycles, one can consider an accelerated nucleon or electron/positron as a particle, which has both charged and neutral states. The acceleration cycle consists of three steps (see Fig. 2). First, a charged particle in relativistic flow is converted into neutral state (point 1). Then, experiencing no influence from the magnetic field, it freely leaves the flow and propagates into ambient medium much further than if it were charged. Second, a transition from neutral to charged state occurs (point 2), which may be the spontaneous neutron decay.

At this moment, the particles in the laboratory frame preserve beaming with the opening angle of $`1/\mathrm{\Gamma }`$, which they had in the neutral state. The initial handicap allows particles to be deflected by an angle $`\theta 2/\mathrm{\Gamma }`$ before the encounter with the relativistic flow. The angular spreading of the particle beam in the laboratory frame means energy gain in the flow comoving frame, which is much larger than $`2`$ and amounts to $`\mathrm{\Gamma }^2`$ in the case of full isotropisation. At the third step, the particles return to the flow (point 3) and their isotropisation in the comoving frame translates into the resulting energy gain in the laboratory frame.

The main parameter, characterizing the efficiency of the converter mechanism, is the optical depth for interactions of accelerated particles (protons/neutrons or electrons/photons) $`\tau =\sigma nD`$, where $`n`$ is the number density of target particles (photons or nucleons) and $`D`$ the accelerator’s size, both measured in the comoving frame, $`\sigma `$ the relevant cross-section.

The optical depth is geometry-dependent. In the case of a continuous outflow or a shock, produced by a central engine and subtending an angle $`>1/\mathrm{\Gamma }`$, one has for photon-induced reactions

$$\tau \frac{\sigma L(\epsilon _{})\mathrm{\Theta }^2}{4\pi Rc\epsilon _{}}.$$
(4)
Here $`L`$ is the apparent luminosity per logarithmic frequency interval at photon energy $`\epsilon _{}=2mc^2\mathrm{\Delta }/(\epsilon \mathrm{\Theta }^2)`$, where the interaction with target photons is the most efficient, $`\epsilon `$ the energy of accelerated particle, $`R`$ the distance from the central engine, $`\mathrm{\Delta }`$ and $`m`$ are the threshold and the mass of the particle for one of the possible cycles.

The beaming angle of target radiation field is $`\mathrm{\Theta }1/\mathrm{\Gamma }`$ in the case where it is produced within the jet or by the shocked gas, and $`\mathrm{\Theta }1`$ for the emission from broad-line regions in AGNs and the radiation scattered in the interstellar medium around GRBs. Intermediate cases, where $`1/\mathrm{\Gamma }<\mathrm{\Theta }<1`$, are also possible. They include, for example, the radiation from inner parts of accretion discs in AGNs or the emission which accompanies the acceleration process (as explained in Sect. V). If the source of target radiation is transient, like a GRB, with duration less than $`R\mathrm{\Theta }^2/c`$, then its luminosity should be replaced by $`Ec/(R\mathrm{\Theta }^2)`$, where $`E=L𝑑t`$.

The most favorable for the converter mechanism conditions exist in AGNs and GRBs, where both ultrarelativistic flows and intense radiation fields are present (see AGNrew  and GRBrew  for reviews). Other objects with relativistic outflows, e.g., stellar-mass microquasars, may also be able to accelerate particles via the converter mechanism. The only two prerequisites are sufficiently high conversion probability (which we specify below) and bulk Lorentz factor $`\mathrm{\Gamma }2`$ in order to compensate the energy losses caused by conversion.

Because of relatively large cross-section, the optical depth for photon-induced reactions is not a bottleneck for the electron cycle, but it could be a limiting factor for the proton cycle. Let us estimate the optical depth for photopionic reactions in three cases. For AGN inner jets we obtain (taking into account only comoving photon fields with $`\mathrm{\Theta }=1/\mathrm{\Gamma }`$)

$$\tau 10^1\left(\frac{L(\epsilon _{})}{10^{45}\text{erg/s}}\right)\left(\frac{\epsilon }{10^{18}\text{eV}}\right)\left(\frac{10}{\mathrm{\Gamma }}\right)^4\left(\frac{10^{15}\text{cm}}{R}\right),$$
(5)
where the apparent luminosity per logarithmic frequency interval $`L(\epsilon _{})`$ depends on $`\epsilon `$. The radiation from AGN broad-line regions creates the optical depth

$$\tau 5\times 10^2\left(\frac{L}{10^{44}\text{erg/s}}\right)\left(\frac{10\text{eV}}{\overline{\epsilon }}\right)\left(\frac{10^{17}\text{cm}}{R}\right),$$
(6)
which does not depend on the particle’s energy: for all particles with energy $`\epsilon >2mc^2\mathrm{\Delta }/\overline{\epsilon }5\times 10^{16}`$ eV the number of target photons is essentially constant because of their relatively narrow spectral distribution. In GRBs, the optical depth due to comoving photons is

$$\tau 3\times 10^3\left(\frac{E(\epsilon _{})}{10^{52}\text{erg}}\right)\left(\frac{\epsilon }{10^{16}\text{eV}}\right)\left(\frac{100}{\mathrm{\Gamma }}\right)^2\left(\frac{10^{16}\text{cm}}{R}\right)^2.$$
(7)
So, the conversion probability $`p_{\mathrm{cn}}=[1\mathrm{exp}(p_c\tau )]`$ in the proton cycle is usually, though not always, much smaller than unity.

It should be noted, that effectively the probability of conversion of a neutron into a proton always exceeds

$$p_{\mathrm{cn}}^{(\mathrm{min})}=\frac{Rm_Nc}{t_n\epsilon }3\times 10^2\left(\frac{10^{15}\text{eV}}{\epsilon }\right)\left(\frac{R}{10^{18}\text{cm}}\right)$$
(8)
because of the neutron decay. Here $`t_n900`$ s is the lifetime of free neutron. The spontaneous decay of free neutrons is important at small energies, especially during the first acceleration cycle, while at large energies the photon-induced conversion is more efficient.

The expressions (5), (6) and (7) cover all physically different situations. For example, one may use Eq. (6) to estimate the optical depth in microquasars, where the target photons are produced by a hot corona having the size $`R10^8`$ cm and the luminosity $`L10^{37}`$ erg/s. The result is $`\tau 0.1`$ at target-photon energy of $`1`$ keV, so that the proton cycle may operate in microquasars at energies $`3\times 10^{14}`$ eV. It is possible, therefore, that some contribution to the galactic CRs around the knee comes from microquasars. In any case, we see no problem in realization of the electron cycle in such objects.

## III Energy gain

Assuming that momenta of particles are isotropized in the comoving frame upon their encounter with the relativistic flow, one gets the average energy gain per cycle

$$g\frac{\left(\mathrm{\Gamma }\theta \right)^2}{2},$$
(9)
where $`\theta `$ is the angle between the particle’s momentum and the flow’s velocity (deflection angle) at the moment of encounter. If the deflection angle is small, $`\theta 1`$, then it grows linearly with distance $`\mathrm{}`$ (travelled by the particle in the charged state)

$$\theta =\theta _0+\mathrm{}/r_g$$
(10)
in uniform magnetic field, while in chaotic field it behaves as

$$\theta =\sqrt{\theta _0^2+\mathrm{}\mathrm{}_c/r_g^2}$$
(11)
on average. Here $`\theta _0`$ is the initial deflection angle, measured at the time of conversion from neutral to charged state, $`\mathrm{}_c<\mathrm{}`$ the turbulence scale of the magnetic field, and $`r_g`$ the gyroradius of the particle, calculated as if the field were uniform. Note, that the turbulence scale is implicitly defined by Eq. (11) and, therefore, may depend on the particle’s energy.

In the case of acceleration at the shock front, the shock catches up to the particle when the displacement of the particle along the shock normal, after it crossed the shock, becomes equal to the distance traveled by the shock front, i.e.,

$$\frac{\sqrt{\mathrm{\Gamma }^21}}{\mathrm{\Gamma }}(\mathrm{}_0+\mathrm{})=\mathrm{}_0\mathrm{cos}\theta _0+_0^{\mathrm{}}\mathrm{cos}\theta d\mathrm{}^{},$$
(12)
where $`\mathrm{}_0`$ is the distance travelled by the particle in the neutral state after it left the relativistic flow, and we neglected the difference between the particle’s velocity and the velocity of light.

In the uniform magnetic field, the deflection angle is

$$\theta \left\{\frac{3\mathrm{}_0}{r_g\mathrm{\Gamma }^2}\left(1\mathrm{\Gamma }^2\theta _0^2\right)\right\}^{1/3}\left(\frac{3\mathrm{}_0}{r_g\mathrm{\Gamma }^2}\right)^{1/3},$$
(13)
provided $`r_g/\mathrm{\Gamma }\mathrm{}_0r_g\mathrm{\Gamma }^2`$. For smaller initial displacement, one arrives to the result of standard theory $`\theta 1/\mathrm{\Gamma }`$, whereas for larger displacement the deflection angle is $`\theta 1`$. In the rightmost part of Eq. (13) we assume that $`\left(1\mathrm{\Gamma }^2\theta _0^2\right)1`$ to simplify the algebra. This corresponds to ignoring particles, which escape propagating (in the comoving frame) nearly parallel to the shock plane. They are neither numerous nor energetically important.

If one substitutes $`\mathrm{}_0`$ by $`R`$ (or by $`\mathrm{\Gamma }D`$ in the case of a small blob with the size $`D<R/\mathrm{\Gamma }`$, ejected by a central engine), then the applicability limits give two critical energies (both measured at the end of acceleration cycle – point 3 in Fig. 2):

$$\epsilon _1=\mathrm{\Gamma }eBR\text{and}\epsilon _2=eBR,$$
(14)
where $`e`$ is the charge of the accelerated particle and $`B`$ the magnetic field strength. The acceleration proceeds with the maximum energy gain of $`\mathrm{\Gamma }^2`$ up to the energy $`\epsilon _2`$, whereas above the energy $`\epsilon _1`$ there is no advantage over the standard mechanism in the energy gain. At the same time, $`\epsilon _1`$ is equal to the maximum energy, achievable in the standard mechanism.

In the chaotic magnetic field

$$\theta \left\{\frac{2\mathrm{}_c\mathrm{}_0}{r_g^2\mathrm{\Gamma }^2}\left(1\mathrm{\Gamma }^2\theta _0^2\right)\right\}^{1/4}\left(\frac{2\mathrm{}_c\mathrm{}_0}{r_g^2\mathrm{\Gamma }^2}\right)^{1/4}$$
(15)
for $`r_g^2/\mathrm{\Gamma }^2\mathrm{}_c\mathrm{}_0\mathrm{\Gamma }^2r_g^2`$, and the critical energies are equal,

$$\epsilon _1=\epsilon _2=\mathrm{\Gamma }eB\sqrt{R\mathrm{}_c}.$$
(16)
In the case of shear flow, the catch-up condition reads

$$\mathrm{}_0\mathrm{sin}\theta _0=_0^{\mathrm{}}\mathrm{sin}\theta d\mathrm{}^{}.$$
(17)
The deflection angle is

$$\theta \left(2\frac{\mathrm{}_0\theta _0}{r_g}\right)^{1/2}\left(\frac{2\mathrm{}_0}{r_g\mathrm{\Gamma }}\right)^{1/2}$$
(18)
for uniform field ($`r_g/\mathrm{\Gamma }\mathrm{}_0r_g\mathrm{\Gamma }`$) and

$$\theta \left(\frac{3}{2}\frac{\mathrm{}_0\mathrm{}_c\theta _0}{r_g^2}\right)^{1/3}\left(\frac{3}{2}\frac{\mathrm{}_0\mathrm{}_c}{r_g^2\mathrm{\Gamma }}\right)^{1/3}$$
(19)
for chaotic field ($`r_g^2/\mathrm{\Gamma }^2\mathrm{}_0\mathrm{}_cr_g^2\mathrm{\Gamma }`$).

There is a subtlety in realization of the converter mechanism if acceleration takes place at a shear flow boundary. If the charged particle reappears at a distance greater than $`r_g`$ (uniform magnetic field) or $`r_g^2/\mathrm{}_c`$ (chaotic field) from the flow boundary, then it should drift or diffuse back. It takes a time of the order of $`t_d=R^2/(\mathrm{\Gamma }r_gc)`$, in the least favourable case of quasi-uniform magnetic field with characteristic spatial scale of $`R`$, or $`t_d=R^2\mathrm{}_c/(\mathrm{\Gamma }^2r_g^2c)`$, in chaotic magnetic field. In both cases we assumed $`\mathrm{}_0\mathrm{sin}\theta _0=R/\mathrm{\Gamma }`$. The solution of diffusion problem with a sink (the shear flow) in the half-space shows, that all particles eventually return to the flow, provided it persists for sufficiently long time. In reality, some of them are lost because the shear flow does not occupy the entire half-space, but the losses are negligible if the spatial extent of the shear flow boundary is much larger than $`\mathrm{}_0\mathrm{sin}\theta _0`$. This condition, for example, is satisfied for a conical jet with the opening angle $`1/\mathrm{\Gamma }`$.

The synchrotron losses during the time $`t_d`$ might be a more serious problem. One can neglect them if the energy loss rate $`\dot{\epsilon }=4/9(e^2/mc^2)^2B^2(\epsilon /mc^2)^2c`$ multiplied by $`t_d`$ is less than the energy of accelerated particle. In the case of quasi-uniform magnetic field, this condition turns out to be energy-independent and, after simple algebra, one gets

$$R\frac{4}{9}\frac{1}{\mathrm{\Gamma }^4}\frac{e^2}{mc^2}\left(\frac{\epsilon _2}{mc^2}\right)^3,\text{i.e.,}R\frac{R_{\mathrm{opt}}}{\mathrm{\Gamma }^2}.$$
(20)
Here $`R_{\mathrm{opt}}`$ is the optimal size of an electromagnetic accelerator prd , which is defined to minimize the amount of energy contained in the electric and magnetic fields. Let us recall, however, that the optimal size is originally defined for the accelerator itself (in this case – for the interior of the shear flow), while Eq. (20) has to do with the external magnetic field.

In a turbulent magnetic field, the synchrotron losses are the most significant for the least energetic particles, i.e., those having $`r_g\mathrm{}_c`$. The losses can be neglected if

$$R\frac{R_{\mathrm{opt}}}{\mathrm{\Gamma }^3},$$
(21)
where the optimal size is the function of critical energy $`\epsilon _2`$, which is defined as if the magnetic field were uniform.

In order to be in agreement with the observed diffuse gamma-ray background, a typical source of the highest-energy CRs must satisfy the condition $`R0.2R_{\mathrm{opt}}`$ prd , so that the synchrotron losses during diffusion do not lead to any further resrictions. The aforementioned statistical arguments are not applicable to individual abberrant accelerators, for which Eq. (20) or Eq. (21) should be considered as an additional limit.

The critical energies in the case of acceleration by a shear flow are

$$\epsilon _1=\epsilon _2=\mathrm{\Gamma }eBR$$
(22)
for uniform magnetic field and

$$\epsilon _1=\mathrm{\Gamma }eB\sqrt{R\mathrm{}_c}\text{and}\epsilon _2=\mathrm{\Gamma }^{\frac{3}{2}}eB\sqrt{R\mathrm{}_c}$$
(23)
for chaotic field. It is interesting to note, that in the latter case $`\epsilon _2>\epsilon _1`$, i.e., the maximum energy achievable in the converter mechanism appears to be $`\mathrm{\Gamma }^{1/2}`$ times larger than that in the standard mechanism. However, a particle can only attain this energy if it enters acceleration cycle with the energy $`\epsilon \epsilon _2/\mathrm{\Gamma }^2`$. A particle, entering the cycle with larger energy, will end up with smaller energy, in contrast to the other three cases.

Of course, every acceleration mechanism must obey the fundamental restrictions and constraints of classical electrodynamics prd . Among them is the Hillas criterion Gr65 ; H84  generalized for relativistic bulk flows, which is simply equivalent to the condition $`\epsilon <\epsilon _1`$, where $`\epsilon _1`$ is the first critical energy from Eqs. (14) and (22). This is not a chance coincidence, but rather the consequence of the implied – uniform – magnetic field configuration. In a turbulent magnetic field, the maximum attainable energy is always lower than the fundamental limits, but the converter mechanism is less affected by the turbulence. Indeed, the magnetic field should be considered chaotic if its turbulence scale is less than the gyroradius of a particle at the beginning of the acceleration cycle. It means $`\mathrm{}_c<r_g(\epsilon _1)`$ for the standard mechanism, but only $`\mathrm{}_c<r_g(\epsilon _2/\mathrm{\Gamma }^2)`$ for the converter mechanism. Thus, a particle accelerated in the standard way feels turbulence starting from larger scale, and a situation is possible, where the magnetic field should be treated as chaotic in the standard mechanism while being essentially uniform for the converter mechanism. The inequality $`\epsilon _2>\epsilon _1`$ from Eq. (23) is, in fact, the consequence of slower decrease of critical energy with decreasing $`\mathrm{}_c`$ for the converter mechanism.

So far, the analysis was limited to the case of negligible synchrotron losses and $`\tau 1`$. The larger is the probability of conversion, the easier is acceleration in the converter mechanism. But if $`\tau `$ approaches unity, then the critical energies decrease: they are still given by Eqs. (14), (16), (22), and (23), but the value of $`R`$ should be replaced by the (energy-dependent) mean free path of the accelerated particle. Analogously, the synchrotron losses also decrease the critical energies; in this case $`R`$ should be replaced by the radiation length.

In both cases the energy share of the accompanying radiation is non-negligible or even dominant, although the status of conversion losses is qualitatively different for the proton and electron cycles (see Sect. V).

## IV Resulting particle distribution

Let us consider a monoenergetic ($`\epsilon =\epsilon _0`$) beam of particles in the neutral state, escaping from the relativistic flow. If the probability of conversion per unit length, $`\lambda `$, is constant (i.e., we neglect the exponential decrease in the number of particles, assuming $`\tau 1`$), then the distribution of charged particles over initial displacement $`dN/d\mathrm{}_0=\lambda `$ translates into

$$\frac{dN}{d\epsilon }\frac{dN}{d\mathrm{}_0}\frac{d\mathrm{}_0}{dg}\lambda \left(\theta \frac{d\theta }{d\mathrm{}_0}\right)^1,$$
(24)
where the particle’s energy is simply proportional to the energy gain, $`\epsilon =g\epsilon _0`$. It is a power-law distribution, but all the particles with displacement larger than the applicability limits of Eqs. (13), (15), (18), and (19) have the constant energy gain of $`\mathrm{\Gamma }^2`$, so that in general a delta-function is added to the resulting particle distribution at its high-energy end. In their power-law parts, the distributions emerging in uniform and chaotic field, respectively, are

$$\frac{dN}{d\epsilon }\epsilon ^{1/2}\text{and}\frac{dN}{d\epsilon }\epsilon $$
(25)
for acceleration at the shock front, and

$$\frac{dN}{d\epsilon }=\mathrm{const}\text{and}\frac{dN}{d\epsilon }\epsilon ^{1/2}$$
(26)
for acceleration at the boundary of shear flow. These spectra are extremely hard. In practice, they can be considered delta-functions, so that the resulting distribution is defined mainly by the spectrum of injected particles.

After a few cycles, a seesaw-shaped spectrum (Fig. 3) is formed from the initial monoenergetic distribution. The number of particles emerging from a cycle, $`N^{}`$, is related to the number of particles entering it as $`N^{}=kN`$, where

$$k=p_{\mathrm{esc}}\underset{i}{}p_{cn}^{(i)}$$
(27)
is the overall probability of conversion, which normally occurs twice per cycle (points 1 and 2 in Fig. 2), multiplied by the probability of escape from the relativistic flow after the first conversion $`p_{\mathrm{esc}}`$, which is the relative number of particles moving upstream with respect to the shock front or the shear flow boundary. Assuming that the particle distribution is isotropic just before the conversion into neutral state, one obtains $`p_{\mathrm{esc}}=1/3`$ for a strong relativistic shock, where the shocked fluid moves at the speed $`c/3`$ away from the shock plane, and $`p_{\mathrm{esc}}=1/2`$ for a shear flow, where the fluid velocity is parallel to the boundary. Within one acceleration cycle, the probabilities $`p_{cn}^{(1)}`$ and $`p_{cn}^{(2)}`$ are of the same order, except for the case where $`p_{cn}^{(2)}`$ is limited by Eq. (8).

There is a simple case, where $`k`$ remains constant in the course of acceleration. This happens in both electron and proton cycles if the spectrum of target photons is a power-law with spectral index $`q=0`$ or (for the proton cycle only) in the case of a narrow-band target-photon field. Then the envelope of the particle distribution shown in Fig. 3 is a power-law $`dN/d\epsilon \epsilon ^\alpha `$, and the index $`\alpha `$ can be obtained in the following way.

Since the particle energies are related as $`\epsilon ^{}=g\epsilon `$, we find that

$$\left(\frac{dN}{d\epsilon }\right)^{}=\frac{k}{g}\frac{dN}{d\epsilon },$$
(28)
where primed is the distribution emerging from the cycle. Taking logarithm of both parts, we obtain

$$\alpha \mathrm{ln}g=\mathrm{ln}\left(\frac{k}{g}\right)\alpha =1\frac{\mathrm{ln}k}{\mathrm{ln}g}.$$
(29)
If the probability of passing through the acceleration cycle $`k`$ is larger than the inverse energy gain $`1/g`$, then the spectral index is $`\alpha <2`$, i.e., most of the energy content is at the high-energy end of the distribution. Both AGNs and GRBs can satisfy this requirement (see Eqs. 5, 6 and 7) and, therefore, can be efficient producers of the highest-energy cosmic rays.

In the energy range, where the efficiency condition for the converter mechanism is satisfied, this mechanism is the dominant source of accelerated particles – just because it provides a spectrum harder than the one resulting from the standard mechanism. At the same time, there is no actual threshold in the conversion probability: the converter mechanism can function even at $`p_{cn}0`$, but the number of the accelerated particles in this case is extremely depleted at high energies. Quantitative studies of the emerging CR spectrum must take into account the dependence of the conversion probability on the energy of accelerated particles. Hence, the precise solution can be obtained only in the self-consistent approach, which includes the effect of the accompanying emission, discussed in the following section.

## V The accompanying emission

Acceleration of particles via the converter mechanism is inevitably accompanied by gamma-ray and – in the case of nucleon conversion – by neutrino emission. The conversion losses in the electron cycle, which are the result of non-optimal inverse Compton scatterings, do not lead to irreversible energy drain. Instead, all the energy remains in photons or electrons (positrons) and, as they both can in principle participate in the acceleration cycle, may be used to inject new particles for acceleration. Thus, even in the case $`\tau 1`$ the influence of the conversion losses is reduced to the decrease of maximum attainable energy and the increase in the level of synchrotron emission by rising the number of particles involved in the acceleration.

On the contrary, the conversion losses in the proton cycle act as a true energy sink: the energy spent for the pion production never comes back to the accelerated nucleons. Neutrinos, $`e^{}e^+`$-pairs and gamma-rays, which are the decay products of pions, are copiously produced as byproducts of the proton acceleration cycle. Neutrinos carry away about one-half of the energy of accelerated nucleons and, because they freely escape from the acceleration site, their spectrum copies that of the nucleons, but scaled down in energies by about an order of magnitude. The gamma-rays and pairs, which carry approximately the same energy as neutrinos, are reprocessed through the electromagnetic cascade: the photons are absorbed in two-photon pair production process, electrons and positrons cool in the magnetic field, producing synchrotron radiation, i.e., another generation of photons, etc.

The hard spectrum of primary electrons and positrons means that they form a standard cooling distribution, $`dN_e/d\gamma _e\gamma _e^2`$, and their synchrotron spectrum is a power-law with (photon) index $`3/2`$. The synchrotron photons spawn another generation of pairs, which cool to form a distribution with index $`5/2`$. The cascade comprises several steps like this, so that the photon spectrum eventually converges to the power-law with index $`2`$. However, at low energies the acceleration site is transparent for the photons, and therefore pairs are not injected below some energy, which leads to a break in the spectrum; below the break the photon spectrum preserves its original spectral index $`3/2`$. The location of the break can be found in a self-consistent way:

$$\epsilon _{\mathrm{br}}0.5\frac{\mathrm{}eB}{\mathrm{\Gamma }m_ec}\left(\frac{\epsilon _{\mathrm{tr}}}{2m_ec^2}\right)^2,$$
(30)
where $`\mathrm{}`$ is the Planck constant. The threshold photon energy $`\epsilon _{\mathrm{tr}}`$ is defined so that the optical depth for two-photon pair production,

$$\tau _{\gamma \gamma }(\epsilon _{\mathrm{tr}})=\frac{\sigma _{\gamma \gamma }}{\sigma _\pi }\tau _\pi \frac{\epsilon _{\mathrm{br}}\epsilon _{\mathrm{tr}}}{2\mathrm{\Gamma }^2m_e^2c^4}\left(\frac{\epsilon _{}}{\epsilon _{\mathrm{br}}}\right)^{1/2},$$
(31)
is equal to unity, where $`\sigma _{\gamma \gamma }`$ is the cross-section of two-photon pair production and $`\tau _\pi `$ the optical depth for photopionic reactions at photon energy $`\epsilon _{}`$. Straightforward calculations yield

$$\epsilon _{\mathrm{br}}0.5\frac{\sigma _\pi }{\sigma _{\gamma \gamma }\tau _\pi }\sqrt{\frac{\mathrm{}eB}{m_ec}\frac{\mathrm{\Gamma }\epsilon _Nm_e^2c^2}{\mathrm{\Delta }_pm_N}}.$$
(32)
Substituting $`B=10^3`$ G, $`\epsilon _N=10^{19}`$ eV (the energy of accelerated nucleons), $`\mathrm{\Gamma }=1000`$, and $`\tau _\pi =10^2`$, – parameters, reasonable for GRBs, – we obtain the break energy $`\epsilon _{br}100`$ keV, which is similar to the really observed in GRBs spectral features. Within this picture, however, there is no simple way to explain the observed relative stability of $`\epsilon _{\mathrm{br}}`$ whereas the parameters entering Eq. (32), especially $`\tau _\pi `$, may vary by orders of magnitude.

Unlike the conversion losses, the synchrotron emission is not tightly related to the acceleration process as such. It may have a negligible effect, especially for protons, but also may be the main energy loss channel. The detailed analysis of the properties of accompanying synchrotron radiation is beyond the scope of this paper, but we point out two distinctive features.

First, the maximum energy of synchrotron photons for the converter mechanism is $`\mathrm{\Gamma }^2`$ times larger than for the standard one. The existence of such an energy limit is easy to see for the standard mechanism (following the arguments of GFR ). The acceleration cycle in this case lasts $`r_g/c`$ and the energy increment is $`\epsilon `$, which gives the acceleration rate $`\dot{\epsilon }\epsilon c/r_g`$. The maximum admissible rate of synchrotron losses is just the same, so that the particle’s energy is limited by the following inequality:

$$\frac{4}{9}\left(\frac{e^2}{mc^2}\right)^2B^2\left(\frac{\epsilon }{mc^2}\right)^2<\frac{\epsilon }{r_g}.$$
(33)
Then, simple calculation yields the maximum energy of synchrotron photons (ignoring relativistic dipole radiation caused by small-scale inhomogeneities of magnetic field):

$$\epsilon _{\mathrm{sy}}^{(\mathrm{max})}0.5\frac{\mathrm{}eB}{mc}\left(\frac{\epsilon }{mc^2}\right)^2\frac{\mathrm{}c}{e^2}mc^2137mc^2.$$
(34)
Doppler boosting gives additional factor $`\mathrm{\Gamma }`$, and in a turbulent magnetic field with the spatial scale of the turbulence less than $`\mathrm{}_c^{\mathrm{cr}}=m^2c^4/\sqrt{e^5B^3}`$ another factor $`(\mathrm{}_c/\mathrm{}_c^{\mathrm{cr}})^{2/3}`$ applies.

The same reasoning is valid for the converter mechanism, in which the cycle duration is $`r_g/c`$ and the energy increment is $`\mathrm{\Gamma }^2\epsilon `$. Consequently, the analog of Eq. (33) gives $`\mathrm{\Gamma }^2`$ times larger limit for the energy of accelerated particle, which translates into factor $`\mathrm{\Gamma }^4`$ in the expression for the energy of synchrotron photons. When the accelerated particle enters the relativistic flow being close to the limiting energy, the synchrotron emission is so efficient, that the particle loses almost all its energy before it is deflected by an angle $`1/\mathrm{\Gamma }`$. Thus, the resulting synchrotron emission is beamed backwards in the flow comoving frame. In the laboratory frame it appears redshifted by the factor $`\mathrm{\Gamma }`$, in contrast to the standard mechanism, in which the synchrotron emission is blueshifted in the laboratory frame by the same factor $`\mathrm{\Gamma }`$. Thus, for an observer resting in the laboratory frame, the maximum energy of synchrotron photons accompanying the converter-acceleration is by the factor of $`\mathrm{\Gamma }^2`$ larger (recall the difference of $`\mathrm{\Gamma }^4`$ in the comoving frame), than in the standard mechanism. Moreover, this highest-energy synchrotron radiation is quasi-isotropic in the laboratory frame, which is another distinctive feature of the converter mechanism. The latter phenomenon has a nature similar to the effect of beam-pattern broadening for the inverse Compton radiation of electrons in front of relativistic shock rtrd . Generally speaking, the converter mechanism makes neutral beams of all kinds (photon, neutrino and neutron beams) broader than $`1/\mathrm{\Gamma }`$, so that they can be seen even if the jet that produced them is not pointing towards the observer.

The accompanying electromagnetic emission of any origin can itself provide photons for the conversion reactions, which may give rise to radiative instabilities analogous the instability driven by $`p+\gamma p+e^{}+e^+`$ process, discussed in BH .

## VI Discussion

The key parameter of the converter mechanism – the probability of conversion – varies from one cycle to another or even within cycles. It increases in the course of acceleration (because a more energetic accelerated particle interacts with less energetic and more abundant target particles) and with increase of the angle between the particle’s trajectory and the velocity of the flow. The first effect is negligible for nucleon-nucleon collisions, while the second is partially or even completely compensated by smaller distance travelled at large angles. There is, however, a general trend of increase of the probability of conversion towards higher energies of accelerated particles; this may even block further acceleration if $`p_{\mathrm{cn}}`$ approaches unity. Nevertheless, with the expansion of the flow the density of target particles drops, leading to the decrease of the probability of conversion, and the acceleration resumes.

It means, that the converter mechanism is capable of self-tuning. The only thing which is required for this mechanism to be efficient is that the probability of passing through the acceleration cycle $`k`$ is larger than $`1/\mathrm{\Gamma }^2`$ somewhere along the flow. It does not matter how large is the probability: in the case $`k1`$ the particles are preaccelerated in the region of high optical thickness and then are further accelerated when (or where) the flow becomes optically thin for them. Of course, the question how much energy is wasted for conversion losses during this intermediate phase is open.

On the other hand, the condition $`k>1/\mathrm{\Gamma }^2`$ places some restrictions on the sources where the proton cycle is readily realizable starting from thermal protons. If a source is known to be bright at hard gamma-rays, one can conclude that the source’s compactness, calculated for the corresponding target photons, is low, so that the probability of proton-to-neutron conversion is $`k\sigma _\pi /\sigma _\gamma 3\times 10^3`$. Thus, the source should have $`\mathrm{\Gamma }20`$ – a condition, satisfied by GRBs and many AGNs. Moreover, this condition is not obligatory for the sources, where the converter mechanism starts from preaccelerated high-energy protons, which require only low-energy target photons.

We expect that the converter mechanism never operates alone and the standard Fermi-type mechanism competes with it. Because the particles spend some time in the neutral state, the converter mechanism has smaller acceleration rate at low energies. But, close to the limiting energy ($`\epsilon _1`$), the durations of acceleration cycle in both mechanisms approach the same value $`R/c`$, so that the average acceleration rates are roughly equal. Except for the case of acceleration in a shear flow with chaotic ambient magnetic field, the mechanisms have just the same absolute, i.e., ignoring radiative losses, energy limit for the accelerated particles. However, in the converter mechanism, this energy is attained in much fewer steps, with potentially much more particles survived. Moreover, the particles – when converted into neutral state – can escape from regions located deep inside the relativistic flow, which further reduces irreversible particle losses in the downstream. So, at the highest energy part of the distribution almost all particles are produced by the converter mechanism, regardless of its performance at low energies.

Now let us consider the phase, when produced CRs leave the accelerator. Within the framework of the converter mechanism, the particles can escape in the form of neutral beam – an easy way, which causes no problem and requires nothing but a sufficiently high conversion probability. In the standard scheme, the escaping particles are charged and inevitably must form an expanding turbulent outflow. This causes considerable adiabatic losses, which can hardly be controlled. Another advantage of the converter mechanism is its greater tolerance for non-uniform magnetic fields, as discussed in Sect. III. Since the magnetic field turbulence has strictly negative effect on the maximum attainable energy, the converter acceleration mechanism may have a larger cut-off energy.

Unlike the standard mechanism, the converter mechanism in many cases does not need any special particle injection. By definition of acceleration cycle, even the particles resting either in the laboratory or in the flow frame serve as injection as soon as they are converted into neutral state. Sometimes however, e.g., in AGNs, the probability of conversion is too small at low energies to provide sufficient injection. Then, the converter mechanism comes to depend on the standard mechanism, which produces preaccelerated particles (at energies $`10^{15}`$ eV in the case of AGNs).

The photopionic processes in AGNs and GRBs leading to copious production of neutral particles have been extensively discussed in the literature (e.g., RM ; AD ), in particular in the context of predictions of fluxes of high energy neutrinos from these objects. The main reasoning was that these objects are postulated to be efficient sources of high-energy CRs and they are surrounded by dense photon fields. Therefore, nucleon-photon interactions may lead to detectable neutrino fluxes. Our point of view is essentially different: if one wants to accelerate CRs in these objects to highest energies, then the converter mechanism is likely the most efficient way to do this, and the neutrino emission appears as a natural and unavoidable byproduct of this mechanism.

Further questions, absent in the test particle problem, arise in the self-consistent approach. In particular, the converter mechanism strongly alters hydrodynamics of relativistic shocks or shear flows. We have shown that they can accelerate nucleons with the efficiency of energy transfer approaching 100%. At such a high efficiency, almost all the available energy is transferred to the end of the particle distribution, and hence the ultra-energetic particles contribute most of the inertia of the relativistic flow (shocked gas). Under such circumstances, the idea of a shock as a discontinuity becomes meaningless, since the gyroradius of dynamically most important particles is comparable to the size of the system.

This limits applicability of the test particle approach, since it is only valid for the particles whose gyroradius is larger than the width of the shock or the shear flow boundary. The similar problem exists also in the standard mechanism, which is, moreover, more sensitive to it. Indeed, in the converter mechanism an accelerated particle may cross the shock or the shear flow boundary while being in the neutral state, hence relaxing the requirements to their sharpness.

The absence of a true shock is an obstacle for the electron cycle, which would otherwise be a very efficient way of direct acceleration of electrons. When the mean free paths of photons and pairs are small compared to the spatial scale of velocity gradient, the direct acceleration of electrons probably cannot rival the efficiency of production of secondary pairs via nucleon acceleration. One of the implications is that at least in some, particularly bright, GRBs the main source of the observed gamma-ray radiation is the pion-induced cascade.

Energetic neutrons (or photons) escape from the relativistic flow and their decay or interaction products can disturb the ambient medium to the extent that it starts to move with an ultrarelativistic speed. An apparent consequence of this is the distortion of beam pattern of the source, while the influence on the limiting energy of the accelerated particles needs further investigation. On one hand, the decrease in the Lorentz factor contrast lowers the limiting energy. On the other hand, enhancement of the ambient magnetic field, either via compression in the shock or by means of various instabilities triggered by neutron decay products in the ambient plasma, has the opposite effect.

## VII Summary

In this paper we suggest and analyze a new acceleration mechanism, which operates via continuous conversion of accelerated particles from charged into neutral state and back. The proposed converter mechanism is efficient for acceleration of both protons and electrons (positrons). It is capable of producing the highest-energy ($`10^{20}`$ eV) cosmic rays in either GRB or AGN environments. With a much lower energy limit, the mechanism can possibly operate in microquasars as well. In the regions of high optical thickness, the converter mechanism is an efficient means of transferring the kinetic energy of bulk relativistic flow to the accompanying radiation, which could explain, for example, the origin of GRB emission. Some peculiarities of the accompanying emission can be a tell-tale sign of the converter mechanism. For example, the production of the highest-energy CRs should be linked with the powerful neutrino emission at a level at least comparable to the power in CRs.

Despite a certain similarity to the standard (diffusive) acceleration mechanism, the converter mechanism violates some of its inherent relations. For example, the maximum particle energy attainable in the converter mechanism is $`\mathrm{\Gamma }^{1/2}`$ times larger than in the standard one, provided acceleration occurs in a shear flow with chaotic ambient magnetic field. Also, the maximum energy of synchrotron photons appears to be $`\mathrm{\Gamma }^2`$ times larger.

Generally speaking, a beam pattern wider than $`1/\mathrm{\Gamma }`$ is characteristic for any type of accompanying emission (synchrotron radiation at highest energies, neutrino emission or photons from pion decay), as well as for the escaping neutrons. This distinctive feature of the converter mechanism opens an interesting possibility for observation of the off-axis blazars and GRBs.

###### Acknowledgements.

We thank John G. Kirk and Michal Ostrowski for useful discussions and the referee for his suggestions. This work was partially supported by the grants 02-02-16236 of the Russian Foundation for Basic Research and 00-15-96674 of the Council for Support of the Leading Scientific Schools, and by programs ”Astronomia” of the Presidium of the Russian Academy of Science and of the Ministry of Industry, Science and Technology.


format: text

text:
# Dynamics of a Complex Quantum Magnet

## Abstract

We have computed the low energy quantum states and low frequency dynamical susceptibility of complex quantum spin systems in the limit of strong interactions, obtaining exact results for system sizes enormously larger than accessible previously. The ground state is a complex superposition of a substantial fraction of all the classical ground states, and yet the dynamical susceptibility exhibits sharp resonances reminiscent of the behavior of single spins. These results show that strongly interacting quantum systems can organize to generate coherent excitations and shed light on recent experiments demonstrating that coherent excitations are present in a disordered spin liquid Ghosh et al. (2002). The dependence of the energy spectra on system size differs qualitatively from that of the energy spectra of random undirected bipartite graphs with similar statistics, implying that strong interactions are giving rise to these unusual spectral properties.

The motivation for understanding the dynamics of quantum systems continues to grow as technology miniaturizes. That quantum and classical dynamics differ fundamentally is underscored by the indications that quantum dynamics can be harnessed in the form of quantum computation to solve problems that are intractable using the processes of classical physics Nielsen and Chuang (2000).

The promise and difficulty of quantum dynamics both arise because the amount of information needed to specify a quantum state grows exponentially with the system size. The number of complex numbers needed to specify a general quantum state of $`N`$ quantum-mechanical two-state spins (or qubits) is $`2^N`$. A classical computer performing $`10^{14}`$ floating point operations per second fas  would take $`10^8`$ years to perform a single operation when $`N=100`$ and $`10^{21}`$ years when $`N=144`$.

Here we study a frustrated spin model that is a quantum generalization of the two-dimensional $`\pm J`$ Edwards-Anderson (E-A) spin glass model Edwards and Anderson (1975), a canonical example of a classical system whose competing interactions give rise to many low-energy states. The essential physical ingredients of the E-A model arise in a wide variety of optimization problems in many fields Mezard et al. (1987): the system cannot satisfy simultaneously all its constraints, and many different configurations are equally effective in minimizing the energy. Though the two-dimensional $`\pm J`$ E-A model is simpler than the three-dimensional version—it is disordered at all nonzero temperatures Palassini and Young (2001) and individual ground states can be found in a time that scales as a polynomial of the system size Bieche et al. (1980)—it has a large ground state degeneracy and a complex energy landscape Palassini and Young (2001); sau ; Santoro et al. (2002). Motivated by recent experiments on $`\mathrm{LiY}_\mathrm{x}\mathrm{Ho}_{1\mathrm{x}}\mathrm{F}_4`$ that demonstrate that quantum tunneling has profound effects on the dynamics of a three-dimensional Ising spin glass with dipolar couplings Brooke et al. (1999, 2001); Ghosh et al. (2002), we study the quantum system obtained by adding a small quantum tunneling term to the two-dimensional E-A model.

We calculate numerically exact eigenstates and dynamical response functions of quantum E-A models with up to $`144`$ spins in the limit of strong interaction strength. At low frequencies the excitation spectrum is remarkably sparse. The number of excitations per unit energy at low energies does not vary significantly with system size, even though the number of eigenvalues grows exponentially with system size while the energy bandwidth is only growing polynomially. The eigenvalues and eigenvectors of each realization are obtained by diagonalizing the adjacency matrix of a graph, and we demonstrate that strong correlations play a vital role by comparing the spin glass excitation spectra to spectra of graphs with similar connectivity statistics but randomly chosen connections. These theoretical results provide a natural framework for understanding recent experiments demonstrating the presence of sharp, saturable resonances in the quantum spin liquid $`\mathrm{LiY}_{0.955}\mathrm{Ho}_{0.045}\mathrm{F}_4`$ Ghosh et al. (2002).

We study two-dimensional systems with periodic boundary conditions in which spin-$`1/2`$ spins interact with nearest neighbors on $`\sqrt{N}\times \sqrt{N}`$ square lattices. The quantum Hamiltonian is

$$H_Q=\underset{ij}{}J_{ij}\sigma _{i,z}\sigma _{j,z}+\mathrm{\Gamma }\underset{i}{}\sigma _{i,x},$$
(1)
where $`\sigma _{i,x}`$ and $`\sigma _{i,z}`$ are Pauli matrices: $`\sigma _{i,x}=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)`$ and $`\sigma _{i,z}=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)`$. The sum $`ij`$ is over all nearest neighbor pairs. A given sample has a fixed realization of bonds in which each bond $`J_{ij}`$ is chosen to be $`J`$ and $`+J`$ with equal probability. This Hamiltonian, which has been studied by many groups tra , is believed to be relevant to the experimental system $`\mathrm{LiY}_\mathrm{x}\mathrm{Ho}_{1\mathrm{x}}\mathrm{F}_4`$ Brooke et al. (1999, 2001); Ghosh et al. (2002), where the coefficient $`\mathrm{\Gamma }`$ is tunable because it is proportional to an applied transverse magnetic field. The $`\mathrm{\Gamma }0`$ limit is the two-dimensional $`\pm J`$ Edwards-Anderson (E-A) model Edwards and Anderson (1975).

As $`\mathrm{\Gamma }/J0`$ the low energy quantum states $`|\psi _n`$ become superpositions of states corresponding to ground state configurations $`|\alpha `$ of the classical model,

$$|\psi _n=\underset{\alpha }{}c_{\alpha n}|\alpha .$$
(2)
The number of ground states of the classical two-dimensional $`\pm J`$ E-A model grows with $`N`$ much more slowly than the number of configurations (though still exponentially) sau ; Landry and Coppersmith (2002); we use an algorithm that we have developed that finds all the classical ground states efficiently Landry and Coppersmith (2002) and standard degenerate perturbation theory Shankar (1980) to compute the low energy quantum states of about $`90\%`$ of $`10\times 10`$ realizations and about $`40\%`$ of $`12\times 12`$ realizations lim . In the limit $`\mathrm{\Gamma }/J0`$, each low energy quantum eigenstate is a superposition of classical ground state configurations related to each other by serial flipping of individual flippable spins, where a flippable spin is one with an equal number of satisfied and unsatisfied bonds Coppersmith (1991). This follows because the matrix element $`\alpha |\left(_i\sigma _{ix}\right)|\beta `$ of the quantum tunneling term between any two classical configurations corresponding to the quantum basis states $`|\alpha `$ and $`|\beta `$ is nonzero only if the two states differ by a single spin flip, and, to lowest order, only states $`|\alpha `$ and $`|\beta `$ with the same energy contribute. The number of flippable spins clearly is no greater than $`N`$, the number of spins in the system, so the matrix characterizing the possible transitions between classical ground states is extremely sparse, and thus well-suited for diagonalization using Lanczos techniques Cullum and Willoughby (1985); Dagotto (1994); Lanczos (1996). We have computed the dynamical magnetic susceptibility by finding low-energy eigenvalues and eigenstates using the sparse matrix ARPACK numerical library Lehoucq et al.  with C++ bindings Shaw . Most of the results shown here are for individual realizations; though there are large sample-to-sample variations in the number of ground states for a given system size Landry and Coppersmith (2002), the qualitative results on which we focus are robust.

The ground state dynamical magnetic susceptibility $`\chi ^{^{\prime \prime }}(\omega )`$ characterizes the response of a system at zero temperature to a magnetic field applied along the $`z`$ axis oscillating at frequency $`\omega `$ sus . The susceptibility consists of sets of Dirac $`\delta `$-function peaks (that in physical systems spread out into a finite width in frequency due to decoherence processes); each peak occurs at a frequency that is $`\mathrm{}`$ times the energy difference between an excited state and the ground state. For the $`\pm J`$ spin glass, as $`\mathrm{\Gamma }/J0`$ the value of $`J`$ affects only the energy zero and the susceptibility at frequencies $`\omega `$ satisfying $`\mathrm{}\omega J`$ depends only on the ratio $`\mathrm{}\omega /\mathrm{\Gamma }`$.

Figure 1 shows the zero temperature dynamic magnetic susceptibility of systems of size $`6\times 6`$, $`8\times 8`$, $`10\times 10`$, and $`12\times 12`$. The density of low energy excitations increases extremely slowly with system size. This result is surprising because the number of energy eigenvalues grows exponentially with the number of spins $`N`$, while these energies all lie within a bandwidth that grows roughly linearly with $`N`$.

To obtain context for these results, we interpret the Hamiltonian matrix for the quantum spin glass as the adjacency matrix of an undirected bipartite graph Chartrand (1985); Cvetkovic et al. (1980); Farkas et al. (2001) in which each classical ground state is a node and edges connect every pair of classical ground states coupled by the quantum term in the Hamiltonian. The graph is bipartite because the edges connect states that differ by a single spin reversal, one of which has an even and the other an odd number of up spins. The spin glass graphs have a modest number of disconnected pieces, called clusters Hartmann (2001). Figure 2 compares the density of energy levels of the largest cluster of a $`10\times 10`$ spin glass realization (with 17040 nodes and 77684 edges) to the density of energy levels of a symmetric bipartite random matrix with 10000 nodes and 50000 edges. The bipartite random matrix has a large energy gap between the ground state and first excited state, and once this gap is exceeded the density of energy levels is much greater than at low energies in the spin glass. The energy level spacing between the excited states of the random matrix is approximately inversely proportional to the number of nodes.

We have compared other properties Albert and Barabasi (2002) of the graphs underlying the quantum spin glass to those of random bipartite graphs. The degree distribution Barabasi and Albert (1999) describing the number of links emanating from the nodes of the spin glass graphs is even narrower than the Poisson distribution of a bipartite random graph with the same mean degree. Figure 3 shows the clustering coefficient $`C`$ Watts and Strogatz (1998), which for bipartite graphs is the probability that two nodes with a common second neighbor are themselves second neighbors Newman et al. (2001). The clustering coefficients of spin glass graphs are significantly larger than those of bipartite random graphs with the same number of nodes and edges Newman et al. (2001), and are close to those of graphs describing $`N`$ noninteracting spins, which have $`2^N`$ nodes, each node with degree $`N`$, and clustering coefficients $`C=4/(N+1)`$clu .

Though some statistical properties of the spin glass graphs are similar to those of graphs for noninteracting quantum spins, the ground state of the quantum spin glass differs significantly from that of noninteracting quantum spins.

Figure 4 (inset) displays the spins in the quantum ground state of a $`12\times 12`$ system that are fluctuating, in that $`|S_{iz}|1`$, where $`S_{iz}`$ is the expectation value of the $`z`$th component of the $`i^{th}`$ spin. The fluctuating spins form connected “bunches” Landry and Coppersmith (2002) with up to 15 spins, so a description in terms of noninteracting spins is not appropriate. The main panel of Figure 4 shows the participation ratios $`P_\alpha `$ of different classical ground states in the quantum ground state, where $`P_\alpha =|c_{\alpha 0}|^2`$, with the $`c_{\alpha 0}`$ defined in Eq. (2). The spin glass participation ratios vary over several orders of magnitude (note the logarithmic abscissa), while for a system of $`N`$ noninteracting spins, the participation ratio is $`1/2^N`$ for all classical configurations. Figure 4 shows that the participation ratios of the adjacency matrix of a bipartite random graph also vary over several orders of magnitude, but with quite different statistics.

We expect the perturbative methods we use to be valid so long as the energy arising from the quantum perturbation ($`<\mathrm{\Gamma }N`$) is smaller than the energy gap between the classical ground state and lowest classical excited states ($`=J`$), so that the procedure is valid only for $`\mathrm{\Gamma }<J/N`$. Exact diagonalizations of very small (3x3 and 4x4) systems are consistent with this expectation.

Our results, particularly the low density of well-defined low-frequency excitations, provide a framework for understanding the generation of coherent, saturable excitations in a complex quantum spin liquid Ghosh et al. (2002). Saturability arises naturally when the density of excitations is low because a frequency that induces a transition from the ground state to an excited state will not be able to induce a second transition from the excited state. When external excitations couple only two states, the system will display coherent oscillations and will have the capability of encoding phase-coherent information Ghosh et al. (2002). An important future goal is to understand how to implement controlled dynamics involving more than two energy levels in this strongly interacting system, as has been achieved in the weakly interacting regime in the context of NMR quantum computation nmr .

We have benefited greatly from discussions with J. Brooke, S. Ghosh, T.F. Rosenbaum, and S.A. Trugman. This work was supported by the MRSEC program of the National Science Foundation under Award No. DMR-9808595 at The University of Chicago and by the National Science Foundation under Award No. DMR-0209630. SNC thanks the Aspen Center for Physics for hospitality during the preparation of this manuscript.


text:
# Non-collinear magnetism in iron at high pressures

## 1 Introduction

Magnetism is known to be important in the phase stability, structure and elastic properties of iron. For example, $`\alpha `$-Fe, the ground state at ambient conditions, would be mechanically unstable if it were not magnetic. Even above the Curie temperature, Tc there are local magnetic moments in $`\alpha `$-Fe. Face-centered cubic iron (fcc or $`\gamma `$-Fe) has incommensurate magnetic correlations which change rapidly with volume, and give rise to the anti-Invar effect (large thermal expansivity) (Mryasov et al., 1992). The magnetic behavior of hcp iron is important for high pressure materials research, and for interpreting high pressure experiments aimed at understanding the Earth’s inner core. Hcp iron is not quenchable to zero pressure, so magnetic studies must be made in situ at high pressures. $`ϵ`$-Fe was long thought to be non-magnetic, due to several Mossbauer experiments that showed no hyperfine splitting in $`ϵ`$-Fe, even down to helium temperatures. However, self-consistent first-principles computations show a magnetic ground stare for $`ϵ`$-Fe which is stable up to about 50 GPa (Steinle-Neumann et al., 1999). Here non-collinear magnetism in $`ϵ`$-Fe is explored using a magnetic tight-binding model fitted to first-principles calculations.

When magnetic moments are collinear, electrons can be considered to be “spin-up” or “spin-down” in a global sense, i.e. there is a global magnetic quantization direction. This means that one can solve for the spin-up and spin-down electrons separately, and then combine the results to compute the total charge density. When moments are oblique to each other, the spin state is said to be “non-collinear.” In that case the problem does not factorize, and one must diagonalize a Hamiltonian of twice the order of the collinear case. A system can be non-collinear either in response to chemical or thermal disorder, or in order to minimize frustration. Both hcp and fcc lattices are frustrated with respect to antiferromagnetism, in that one cannot order these lattices antiferromagnetically so that all neighbors have opposite spins. This is known to lead to spin-waves and non-collinear magnetism in fcc-Fe. Also, a common example of non-collinear spins occurs on heating bcc-Fe above the Curie point, where the spins disorder dynamically. A comprehensive review of the theory of non-collinear magnetism is given in Sandratskii (1998).

## 2 Method

To study non-collinear magnetism in Fe, a first-principles based non-magnetic tight-binding model (Cohen et al., 1994, 1997; Wasserman et al., 1996) is combined with a model for magnetism (Pickett, 1996; Mukherjee and Cohen, 2001). The Hamiltonian (H) of the system is given by,

$$H=H_0+R^+H_{collinear}R,$$
(1)
where $`H_0`$ is the (doubled) non-magnetic Hamiltonian, $`H_{collinear}`$ is the collinear magnetic Hamiltonian and $`R`$ is the rotation matrix for spin directions, as described below. The overlap matrix $`S`$ is unchanged from the non-magnetic non-orthogonal tight-binding model (Cohen et al., 1994, 1997; Wasserman et al., 1996). $`H_0`$ was fit to eigenvalues and total energies from an extensive set of non-magnetic Linearized Augmented Plane Wave (LAPW) results within the Generalized Gradient Approximation (GGA) (Perdew et al., 1996) in the bcc, fcc, and hcp structures as functions of pressure and strain, and has been extensively tested (Cohen et al., 1994, 1997; Wasserman et al., 1996).

The magnetic Hamiltonian is twice the size of the original non-magnetic Hamiltonian, with the upper left block being spin-up, the lower right block is spin down, and the off-diagonal blocks are from coupling between up and down spins that give rise to non-collinear magnetic solutions. $`H_0`$ is doubled from the original non-magnetic Hamiltonian, with identical diagonal blocks, and zeros on the off-diagonal blocks. The only non-zero elements in the collinear magnetic Hamiltonian Hcollinear are the diagonal onsite $`d`$ elements. The latter are given by $`Im_i/2`$ for the spin-up block, and $`+Im_i/2`$ for the spin-down block, where $`m_i`$ is the magnetic moment of the atom $`i`$, and $`I`$ is the Stoner parameter, which controls the strength of the exchange splitting (Pickett, 1996). $`R`$ is the rotation matrix which depends upon the atomic spin direction with respect to a global reference frame (Uhl et al., 1994). The diagonal of $`H_{mag}`$ is spin-up and spin-down, and the off-diagonal bands are the coupling between up and down, which arise from the rotation operation. Similar models have been used previously for incorporating magnetism within a tight-binding approach (You and Heine, 1982; Pickett, 1996; Freyss et al., 1997; Mehl et al., 2001), but the magnetic behavior of $`ϵ`$-Fe has not been addressed, nor has a first-principles tight-binding model previously been developed that includes constraining fields.

The full spin density matrix is computed, and the constraining fields are computed self-consistently so that spin directions as well as spin moments are self-consistent. The spin-density matrix is given by

$$\rho _{\sigma \sigma ^{}}=\underset{ij}{}\psi _{i\sigma }S_{ij}\psi _{j\sigma ^{}}$$
(2)
where $`S_{ij}`$ is the overlap matrix and $`\psi _{i\sigma }`$ are the eigenvectors for atom orbital $`i`$ and spin $`\sigma `$ in the local coordinate system on each atom. In the local coordinate system the spin on each atom is diagonal and can be represented as up or down. The spin-1/2 rotation matrix elements of $`R`$ that rotates from the local coordinate system to the global system with oblique spins is given by

$$R=\left(\begin{array}{cc}R_{}& R_{}\\ R_{}& R_{}\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\frac{\theta }{2}e^{i\frac{\varphi }{2}}& \mathrm{sin}\frac{\theta }{2}e^{i\frac{\varphi }{2}}\\ \mathrm{sin}\frac{\theta }{2}e^{i\frac{\varphi }{2}}& \mathrm{cos}\frac{\theta }{2}e^{i\frac{\varphi }{2}}\end{array}\right)$$
(3)
where $`\theta `$ is the polar angle and $`\varphi `$ is the azimuthal angle. The matrix elements given here are used to build the full rotation matrix; the only non-zero elements are diagonal in atom orbital indices. Without spin-orbit coupling, there is no coupling to absolute directions relative to the crystalline lattice, and only the relative angles between spins are important. Spin-orbit interactions are not included here; their energetic importance is insignificant compared with the exchange energy changes for iron.

The magnetic moments are given by the spin density on each atom summed over the orbitals. The expressions are

$`m`$ $`=`$ $`\sqrt{(\rho _{}\rho _{})^2+4\rho _{}\rho _{}^{}}`$ (4)
$`\theta `$ $`=`$ $`\mathrm{arccos}\left({\displaystyle \frac{\rho _{}\rho _{}}{m}}\right)`$ (5)
$`\varphi `$ $`=`$ $`\mathrm{}\mathrm{log}\rho _{}.`$ (6)
The magnetic moments are determined self-consistently. One starts with an initial guess, which gives the exchange splittings, $`\pm Im_i/2`$, in the Hamiltonian, find $`H`$, and diagonalize the generalized eigenproblem. States are occupied up to the Fermi level, and the output moment on each atom is found from the weighted eigenvectors and the expressions above. The process is repeated until self-consistency is obtained.

In general the output moment direction, as well as magnitude, will differ from the input. In this way a self-consistent process can be used to find the moment magnitudes and directions that are locally stable. In some cases one wants to compute the energy of a given magnetic structure. In order to guarantee that the output magnetic structure is the desired structure, it is necessary to apply magnetic fields that force the output moments to form the desired structure. Such a procedure has been implemented using a self-consistent procedure with the option to constrain the moments in magnitude and direction (constraining fields contain longitudinal components), or to allow the magnitude to adjust self-consistently, but force the direction to remain as desired (constraining fields are transverse). The procedure used is similar to, but different from, the procedure outlined in Ujfalussy et al. (1999).

The total energy is given by

$$E_{total}=E_b+Im_i^2/4+\stackrel{}{b}_i\stackrel{}{m}_i$$
(7)
where $`E_b`$ is the band energy (sum of the eigenvalues). There are no additional potential terms; the total energy is obtained from the band structure and moments. The second term in Eq. 7 corrects for double counting of the exchange interaction (Pickett, 1996), and the third term corrects for double counting of the interaction with the applied staggered fields $`\stackrel{}{b}_i`$.

The model differs from the conventional Stoner model (Stoner, 1938), in which a ferromagnetic instability is predicted by the inequality $`IN(0)>1`$, where $`N(0)`$ is the non-magnetic density of states at the Fermi level, and the extended Stoner model (Krasko, 1987), in which $`N`$ is replaced by the effective density of states $`\stackrel{~}{N}(M)=M/\delta ϵ`$, where $`\delta ϵ`$ is the exchange splitting, in that our model allows for different hybridization depending on magnetic state, accounting thereby to the actual magnetic structure, i.e. whether the system will be ferromagnetic, antiferromagnetic, or non-collinear.

## 3 Results and Discussion

### 3.1 bcc

The model is in good agreement with previous self-consistent calculations for $`\alpha `$-Fe, especially when a volume dependent $`I`$ is used (Mukherjee and Cohen, 2001). The lowest energy state for $`\alpha `$\- Fe for the model at ambient and at higher pressures is ferromagnetic, in agreement with the first principles LAPW self-consistent calculations (Stixrude et al., 1994). Moreover, even with a single value of volume independent $`I`$ the model calculation describes well the pressure dependence of the magnetic moment and the magnetization energy. As pressure is increased the magnetization energy decreases smoothly (Fig. 1). The value of $`I`$ needed for quantitative agreement is almost constant (around 1 eV) for larger atomic volumes or lower pressures but increases for higher pressures. This is consistent with the increase in the exchange interaction at higher pressures noted previously (Asada and Terakura, 1992). Even though the exchange interaction increases with pressure, the non-magnetic density of states near the Fermi energy decreases, giving rise to a net reduction in the magnetic moment and magnetization energy, and finally the loss of magnetism at around atomic volume of 40 bohr3.

The Stoner parameter $`I`$ was adjusted at each volume to give agreement with self-consistent full-potential Linearized Augmented Plane Wave (LAPW) magnetization energies (Steinle-Neumann et al., 1999) for each volume for bcc. $`I`$ is expected to increase as the electronic density increases with decreasing atomic volume, similar to what has been seen before for $`\gamma `$-Fe (Krasko, 1987) and $`\alpha `$-Fe (Fig. 1). Fitting the results to a polynomial valid between $`V=`$50 bohr3 and 90 bohr3 gives $`I(eV)=3.41260.08583V+0.00103V^24.166610^6V^3`$ .

### 3.2 fcc

$`\gamma `$-Fe shows a great richness in non-collinear magnetic structures as volume is varied, and here it is not explored in detail; much study has been done, as reviewed in Sandratskii (1998). The energies and moments for fcc-Fe with spiral spins along (001) and $`\theta =\pi /2`$ were determined using eight-atom supercells, allowing calculations for $`(00q)`$ with $`q=0,\pi /4,\pi /2,3\pi /4`$,and $`\pi `$. A k-point mesh of 12$`\times `$12$`\times `$4 was used, giving 72 k-points in the irreducible wedge for tetragonal symmetry. For $`\gamma `$-Fe, the complex behavior of magnetic moments and the magnetization energy with increasing pressures and varying spiral spin density wave states is qualitatively reproduced by the model (Figs. 2 and 3) compared with self-consistent calculations (Mryasov et al., 1992; Uhl et al., 1994; Sjöstedt and Nordström, 2002). Quantitatively the results are sensitive to the value of $`I`$, and to get good agreement with self-consistent calculations it seems a smaller I, about 0.94 eV, is required than derived from bcc-Fe (0.99-1.02 eV for the volumes considered for fcc). In any case there is some variation in self-consistent results for $`ϵ`$-Fe, due to extreme sensitivity of the magnetic structure to basis set, k-point sampling, etc., due to the small energy scale. Furthermore, the magnetic ground state in $`\gamma `$-Fe is sensitive to the atomic moment approximation (Sjöstedt and Nordström, 2002) in the tight-binding model; that is the moment is really a field that varies with position in the crystal, and is not constant on each atom. Within the atomic moment approximation our results are reasonably consistent with the self-consistent results. More detailed comparisons with the comprehensive non-collinear LAPW calculations of Sjöstedt and Nordström (2002) including k-point convergence tests are called for, but have not yet been done. In Sjöstedt and Nordström (2002) 4000 k-points were used in the full Brillouin zone, compared to our 4608, which seems comparable. It is more difficult to compare with Mryasov et al. (1992) since a real space multiple scattering approach and a muffin-tin potential approximation was used.

### 3.3 hcp

Previous collinear first principles calculations show stability of an antiferromagnetic state, afmII, in hcp iron (Steinle-Neumann et al., 1999). Moreover, the computed equation of state of afmII greatly improves agreement with the experimental equation of state. An hcp lattice is frustrated for antiferromagnetism; it is not possible to have perfect antiferromagnetic order on it. In the afmII structure each atom has eight antiferromagnetically oriented and four ferromagnetically oriented neighbors, maximizing the antiferromagnetic interactions. In the case of a nearest-neighbor (n.n.) Heisenberg model, which has energy $`E=J_1_{n.n.}\stackrel{}{m}_i\stackrel{}{m}_j`$, the energy is independent of the angle between the moments of one antiferromagnetic pair and another (see Fig.4). If non-neigherest neighbor interactions are important, or if the Heisenberg model does not completely describe the energetics, the energy might be further lowered if the antiferromagnetic pairs are oblique or perpendicular to each other (i.e. $`\alpha 0`$).

The collinear afmII structure has 4 atoms per unit cell in space group Pmma (Steinle-Neumann et al., 1999). In order to tile the lattice with the pattern shown in Fig.4 a cell with 8 atoms is obtained with space group A value of 1.6 was used for $`c/a`$ for all of the hcp based calculations here, which is close to the ground state value for the volumes studied here. Tests showed that $`c/a`$ does not vary significantly (¡.005) with magnetic state. Figure 5 shows the energy versus angle $`\alpha `$ for different volumes. In all cases the energy decreases when $`\alpha `$ is varied from zero, indicating that the ground state is non-collinear. The most stable state is with $`\alpha `$=90, which minimizes the local ferromagnetic interactions. Table 1 shows the results of the TB model for hcp Fe with ferromagnetic, admII, and the non-collinear structure with $`\alpha `$=90 (ncl).

The ground state non-collinear total energies of $`ϵ`$ -e were fit to a Vinet equation of state (Vinet et al., 1987) as a function of volume. Table 2 shows the results both using $`K^{}`$=4 and allowing $`K^{}`$ to vary. For non-magnetic Fe, $`K^{}`$=4 is not a bad approximation, but when magnetism is included $`K^{}>`$6. The results show that the TB model is in god agreement with the self-consistent LAPW computations of Steinle-Neumann et al. (1999). Including magnetism lowers the bulk modulus, and including non-collinear magnetism lowers it further. The resulting equation of state fo ncl is in very good agreement with the experiments (Fig. 6,Table 2) (Jephcoat et al., 1986; Mao et al., 1990). When magnetization is not included the disagreement between the experiment and theory is 75$`\%`$ and 9$`\%`$ for the bulk modulus and the equilibrium volume respectively. Including the afmII structure as the ground state of hcp reduces the disagreement to 25$`\%`$ and 5$`\%`$ for bulk modulus and equilibrium volume (Steinle-Neumann et al., 1999). When the non-collinear spin state energies were used in the equation of state the bulk modulus was within 10$`\%`$ of the experiment.

In spite of theoretical evidence for magnetism in $`ϵ`$-Fe, and the great improvement in the equation of state when magnetism is included, the experimental situation is unclear. Mössbauer spectroscopy shows no evidence of magnetism in hcp Fe (Taylor et al., 1982, 1991). X-ray absorption spectroscopy also has been interpreted to show lack of magnetism (Rueff et al., 1999), but due to the absence of an absolute calibration, and sensitivity of the spectrum to changes in the density of states, the experiments show only that magnetic moments are lower in $`ϵ`$-Fe than $`\alpha `$-Fe, a result consistent with theory. On the other hand, there is some independent evidence of magnetism from Raman spectroscopy, which shows two peaks (Merkel et al., 2000), instead of the one expected in the hcp structure. The frequencies and splitting of these peaks is predicted well from first-principles calculations for magnetic ordered afmII $`ϵ`$-Fe (Steinle-Neumann et al., 2003), and the absence of observed splitting in Mössbauer will be explained in a subsequent paper.

Thanks to B. Fultz, A. Goncharov, R. Hemley, A. Jephcoat, I. Mazin, H. Olijnyk, W. Pickett, and G. Steinle-Neumann, L. Stixrude, V. Struzhkin for helpful discussions. S. Mukherjee performed technical assistance. This work was supported by DOE ASCI/ASAP subcontract B341492 to Caltech DOE W-7405-ENG-48 and the National Science Foundation EAR-998002. Computations were performed on the Cray SV1 at the Geophysical Laboratory, supported by NSF grant EAR-9975753 and by the W. M. Keck Foundation.