text:
# Microcanonical determination of the order parameter critical exponent

## I Introduction

Substantial progress has been achieved in the understanding of discontinuous phase transitions since the problem has been formulated in microcanonical terms. Till now the microcanonical analysis of continuous phase transitions has only been partly successful. On the one hand it is a great accomplishment that the typical features of symmetry breaking, i.e. the abrupt onset of the order parameter as the critical point is approached from above and the diverging susceptibilities, turn up already for finite systems. This is in contrast to the canonical ensemble where singularities appear exclusively in the thermodynamic limit. On the other hand it is an intriguing and disappointing fact that mean field values have been found for all critical exponents and for all finite system sizes when directly analysing the density of states.

We do not contest the earlier analysis where values of the critical exponents of an Ising system were determined from the expansion of the entropy $`S_N(E,M)`$ in the vicinity of the critical point of the finite system, situated at $`E=E_{N,c}`$ and $`M=0`$. Here it is our aim to show that the non classical infinite lattice exponents can be obtained directly from the density of states if the expansion of the finite system entropy $`S_N(E,M)`$ is carried out at the critical point $`E=E_{c,\mathrm{}}`$ and $`M=0`$ of the infinite system. (The critical point of the finite or infinite system is defined as the value of the energy $`E=E_{c,N}`$, where the second derivative $`[^2S_N/M^2]_E`$ of the entropy taken at $`M=0`$ changes its sign.)

For the purpose of determining the density of states $`\mathrm{\Omega }_N(E,M)=\mathrm{exp}S_N(E,M)`$ (in units where $`k_B=1`$) we have developed a highly efficient algorithm which is based upon the Broad Histogram method. It allows a precise and speedy determination of the density of states as a function of one, two or even more parameters should it be necessary. In the course of these calculations it will be demonstrated that very good values for the order parameter critical exponent can be obtained by a microcanonical analysis already for relatively small system sizes.

The paper is organized in the following way. In the next Section we present the numerical method which permits us to compute the density of states of classical spin systems with very high accuracy. Section III is devoted to the analysis of the microcanonically defined spontaneous magnetization and to the determination of the corresponding critical exponent. A short summary concludes the paper.

## II Numerical method

One of the aims of computer simulations in statistical mechanics is the determination of averages of observables $`\widehat{X}`$ of the many particle system under study:

$$\widehat{X}_{\widehat{W}}=\underset{\mu \mathrm{\Gamma }_\mu }{}\widehat{W}(\mu )\widehat{X}(\mu ).$$
(1)
The space $`\mathrm{\Gamma }_\mu `$ contains an enormous number $`\mathrm{\Lambda }_\mu `$ of microstates and the sum in (1) runs over all of them. $`\widehat{W}(\mu )`$ is a probability distribution on that space and $`\widehat{X}(\mu )`$ is the value of $`\widehat{X}`$ in the microstate $`\mu `$.

In all cases of interest $`\widehat{W}(\mu )`$ and $`\widehat{X}(\mu )`$ depend on $`\mu `$ only via a few parameters $`\underset{¯}{P}=(P_a,P_b,\mathrm{})`$, i.e.,

$$\widehat{W}(\mu )=W(\underset{¯}{P}(\mu ))$$
(2)
and

$$\widehat{X}(\mu )=X(\underset{¯}{P}(\mu ))$$
(3)
Examples for $`\underset{¯}{P}`$ might be the interaction energy $`E`$ of a magnetic system and its magnetization $`M`$. An example for $`W(\underset{¯}{P})`$ could be the canonical distribution

$$W(E,M)=\frac{\mathrm{exp}\left[\beta \left(EhM\right)\right]}{\underset{E^{},M^{}}{}\mathrm{exp}\left[\beta \left(E^{}hM^{}\right)\right]}.$$
(4)
The parameters $`\underset{¯}{P}`$ characterize the macrostate of the system and examples for $`X(\underset{¯}{P})`$ might be the parameters themselves, their fluctuations as e.g. $`(E<E>_W)^2`$ or any other function of them.

The unwieldy sum in (1) is greatly reduced if the degeneracies $`\mathrm{\Omega }(\underset{¯}{P})`$ of the macrostates which are characterized by the parameters $`\underset{¯}{P}`$ are known. Then:

$$\widehat{X}_{\widehat{W}}=X_W=\underset{\underset{¯}{P}}{}\mathrm{\Omega }(\underset{¯}{P})W(\underset{¯}{P})X(\underset{¯}{P}).$$
(5)
$`\mathrm{\Omega }(\underset{¯}{P})`$ is also called the density of states DOS. The number of terms in the sum (1) is already gigantesque for a system of only a few hundred particles: $`\mathrm{\Lambda }_\mu `$ is of order $`f^N`$ where $`N`$ is the number of particles in the system. $`f`$ is equal to 2 for an Ising system and for a $`q`$-states Potts model $`f=q`$.

Compared to that, the sum in (5) can be performed with very little effort indeed. Its number of terms is only of order $`N^\nu `$ where $`\nu `$ is the number of parameters. Once $`\mathrm{\Omega }(\underset{¯}{P})`$ is known, the sum in (5) can be performed in the twinkle of an eye by any desktop computer.

Now the difficulty of calculating averages of $`X`$ has been shifted to the determination of the DOS. This does not seem to be any easier, but we have developed a new algorithm to construct $`\mathrm{\Omega }(\underset{¯}{P})`$. Before the algorithm is explained in detail we want to point at yet another advantage which is offered at no extra cost when the density of states $`\mathrm{\Omega }(E,M)`$ has been determined. As $`\mathrm{ln}\mathrm{\Omega }(E,M)=S(E,M)`$ (setting $`k_B=1`$) is the entropy as a function of its natural variables, the thermodynamic relations can be derived directly from $`S(E,M)`$ by differentiation. One thereby obtains the microcanonically defined temperature, specific heat, susceptibilities etc.

This is an extra benefit on top of the possibility of having the mean values of any function $`X(\underset{¯}{P})`$ (e.g. $`X(E,M)`$) and for any desired probability distribution $`W(\underset{¯}{P})`$ at hand.

The algorithm is a variant of the transition variable method. Its starting point is an extremely simple observation: Consider any mechanism which generates a new microstate $`\mu ^{}`$ when it is applied on the microstate $`\mu `$. An example for such a mechanism could be a single spin flip in an Ising system. When the single spin flip mechanism is applied to all the $`N`$ spins of all the $`\mathrm{\Omega }(\underset{¯}{P}_i)`$ microstates which belong to the macrostate which is characterized by $`\underset{¯}{P}_i`$, then $`N\mathrm{\Omega }(\underset{¯}{P}_i)`$ new microstates are created. Let $`V_{ij}`$ be the number of them which belong to the macrostate characterized by $`\underset{¯}{P}_j`$. As any spin flip can be reversed: $`V_{ji}=V_{ij}`$. In other words there are as many connections (by single spin flip) from level $`\underset{¯}{P}_i`$ to level $`\underset{¯}{P}_j`$ as there are connections from level $`\underset{¯}{P}_j`$ to $`\underset{¯}{P}_i`$. Thus if one chooses one microstate of level $`\underset{¯}{P}_i`$ at random, selects one spin at random, and flips it, the probability of arriving at level $`\underset{¯}{P}_j`$ will be $`V_{ij}/(N\mathrm{\Omega }(\underset{¯}{P}_i))`$ and vice versa if we choose one of the spins in one of the states of level $`\underset{¯}{P}_j`$ at random and apply the mechanism to it, then the probability of arriving at level $`\underset{¯}{P}_i`$ is $`V_{ij}/(N\mathrm{\Omega }(\underset{¯}{P}_j))`$.

This simple observation is used to create an algorithm which serves two purposes:

i) It determines the rate of attempts $`T_{ij}=V_{ij}/(N\mathrm{\Omega }(\underset{¯}{P}_i))`$ to go from a starting level $`\underset{¯}{P}_i`$ to another level $`\underset{¯}{P}_j`$ which can be reached by the chosen mechanism.

ii) It leads to an equal probability of visiting the levels irrespective of their degeneracy.

Once the $`T_{ij}`$ have been found it is an easy matter to construct the whole surface $`\mathrm{\Omega }(\underset{¯}{P})`$ from the ratios $`T_{ij}/T_{ji}=\mathrm{\Omega }(\underset{¯}{P}_j)/\mathrm{\Omega }(\underset{¯}{P}_i)`$.

In order to calculate the $`T_{ij}`$ one writes down the total number $`Z_i`$ of applications of the mechanism while in level $`\underset{¯}{P}_i`$ and also the number of applications $`Z_{ij}`$ which would lead to level $`\underset{¯}{P}_j`$ and this in both cases, if the step is executed or not. Then $`T_{ij}=Z_{ij}/Z_i.`$

As a reversible transition mechanism we choose single spin flips. Before the spin is flipped the system is characterized by the parameter values $`\underset{¯}{P}_i`$. A spin is chosen at random. Its flip would lead to a new state characterized by $`\underset{¯}{P}_j`$. We add 1 to $`Z_i`$ and also to $`Z_{ij}`$ whether the flip is accepted or not. After a large number of attemped spin flips $`Z_{ij}/Z_i`$ approches $`T_{ij}`$. The demand ii) is fulfilled, if steps from $`i`$ to $`j`$ and from $`j`$ to $`i`$ are executed with the same frequency. This can be achieved by introducing probabilities for the acceptance of executing a step. This probability of accepting the step from $`i`$ to $`j`$ is equal to 1 if $`T_{ij}<T_{ji}`$ and it is $`T_{ji}/T_{ij}`$ otherwise. This equalizes the probabilities

$$p(\underset{¯}{P}_i,\underset{¯}{P}_i+\mathrm{\Delta }\underset{¯}{P})=p(\underset{¯}{P}_i+\mathrm{\Delta }\underset{¯}{P},\underset{¯}{P}_i)$$
(6)
of going from a level $`\underset{¯}{P}_i`$ to the level $`\underset{¯}{P}_i+\mathrm{\Delta }\underset{¯}{P}`$ and of the reverse transition.

Choosing nonzero, but otherwise arbitrary values for $`Z_i`$ and $`Z_{ij}`$ for all $`i`$ and $`j`$ and starting from any one of the states of the system one rapidly builds up good estimates for the transition rates $`T_{ij}`$. When this point is reached all levels $`\underset{¯}{P}_i`$ are visited with the same frequency and the quality of the transition variables is improved with the same rate over the whole parameter range.

The procedure adopted here differs from other implementations of the transition variable method in several respects: The transition variables are updated at every single spin flip, in return we do not register all the possible transitions which are possible from a given microstate, but only those ones which are attempted by the algorithm, whether they are successful or not. Furthermore the transition probabilities $`T_{ij}`$ are updated during the whole run and the acceptance rates are governed by the most recent value of $`T_{ij}`$ from the beginning of the simulation to its very end.

This method gives us the freedom of restricting the calculations to a chosen number of channels in energy or magnetization direction. Whereas for not too large systems all values of the energy and of the magnetization may be allowed in a single run, for larger system sizes narrow bands in energy direction covering all magnetization channels may be used. In the limit of only one energy channel per band our approach differs from the well-known Q2R method , as for a fixed energy value all the possible magnetizations are visited with the same frequency.

Our algorithm has been subjected to several stringent tests. We have calculated $`\mathrm{\Omega }_N(E)`$ for a 32 x 32 Ising model and we have compared the resulting entropy with the exact data of Beale. $`S_N(E)`$ as well as its first and second derivate match the exact values almost perfectly. For example, $`7\times 10^5`$ Monte Carlo sweeps yield for the inverse microcanonical temperature $`S_N^{}(E)`$ an average error as small as $`0.025\%`$ when compared to the first derivative of the Beale data. In a plot of the specific heat $`c(E)=(S_N^{}(E))^2/S_N^{\prime \prime }(E)`$ versus temperature $`T(E)`$, where $`1/T(E)=\beta ^{}(E)=S_N^{}(E)`$, one cannot distinguish the simultated data from the exact result. For the same system we have also calculated $`\widehat{\mathrm{\Omega }}_N(E,M)`$ from which we again find $`\mathrm{\Omega }_N(E)`$, now by summation over M. In addition we find $`\stackrel{~}{\mathrm{\Omega }}_N(M)`$ by summation over $`E`$ and there is exellent agreement with the exact result which is $`N!/[\frac{N+M}{2}!\frac{NM}{2}!]`$.

The algorithm compares very favorably to other methods proposed for the computation of the density of states. Its efficiency for the 32 x 32 Ising model is slightly better than that of the multiple-range random walk algorithm of Wang and Landau where an average error of $`0.035\%`$ was obtained for $`7\times 10^5`$ Monte Carlo updates. It must be stressed that our method is completely different to the Wang/Landau method where a rough estimator of the DOS is built up rapidly due to a multiplicative update of $`\mathrm{\Omega }(E)`$ every time the energy level $`E`$ is visited. The error introduced by this multiplicative update has to be eliminated before obtaining the DOS with high precision . In our present approach, no such error is introduced during the course of the simulation.

As we will show in the following, the new method enables us to compute the density of states with a high accuracy even for large Ising systems. Furthermore, preliminary results for complex systems with a rough energy landscape are very promising .

As a first application, we compute in the next Section the microcanonical order parameter in the two- and three-dimensional Ising models. It must be noted that data of extremely high accuracy, rarely needed in a canonical study, are mandatory for a reliable determination of microcanonically defined quantities.

## III Microcanonically defined spontaneous magnetization

In a microcanonical analysis of finite systems with $`N=L^d`$ spins ($`d`$ being the number of space dimensions), the object of interest is the density of states $`\mathrm{\Omega }_N(E,M)`$ as a function of the energy $`E`$ and the magnetization $`M`$. We define the spontaneous magnetization $`M_{sp,N}(E)`$ as the value of $`M`$ where the entropy $`S_N(E,M)=\mathrm{log}\mathrm{\Omega }_N(E,M)`$ at a fixed value of the energy has its maximum with respect to $`M`$. At energies lower than a critical value $`E_{c,N}`$ the entropy exhibits two maxima at $`M=\pm M_{sp,N}`$. On approach to this point $`E_{c,N}`$ the order parameter $`M_{sp,N}`$ vanishes with a square root behaviour in the finite system .

In order to achieve very good statistics we have restricted the calculations to a narrow band which covers only five channels in energy direction, but stretches over all possible values of the magnetization. All single spin flips within this band are allowed. Spin flips which would end up outside the band are rejected. We build up the $`T_{ij}`$ and from them construct the DOS in magnetization direction. A single flip leads from a microstate in the level $`\underset{¯}{P}_i=(E_i,M_i)`$ with entropy $`S_N(E_i,M_i)`$ to a state of a neighbouring level $`\underset{¯}{P}_j=(E_j,M_j)`$ with $`E_j=E_i+4\eta `$ and $`M_j=M_i+2\gamma `$. Here $`\gamma =\pm 1`$, whereas $`\eta =0`$, $`\pm 1`$, $`\pm 2`$ for a two-dimensional square lattice and $`\eta =0`$, $`\pm 1`$, $`\pm 2`$, $`\pm 3`$ for a three-dimensional simple cubic lattice. The entropy at the macro state at $`\underset{¯}{P}_j`$ can then be expressed as:

$`S_N(E_j,M_j)=`$ $`S_N(E_i,M_i)+\beta ^{}(E_i,M_i)\mathrm{\hspace{0.17em}4}\eta `$ (8)
$`+\mu (E_i,M_i)\mathrm{\hspace{0.17em}2}\gamma +\mathrm{}`$
where $`\beta ^{}(E,M)=S_N(E,M)/E`$ is the inverse temperature and $`\mu (E,M)=S_N(E,M)/M`$. Because of the densely spaced energy and magnetization levels (on the intensive scale) higher order terms in (8) may savely be neglected. From ratio of one transition observable pair $`T_{ij}/T_{ji}=\mathrm{\Omega }(\underset{¯}{P}_j)/\mathrm{\Omega }(\underset{¯}{P}_i)=\mathrm{exp}\left[4\eta \beta ^{}+2\gamma \mu \right]`$ follows one linear combination of $`\beta ^{}`$ and $`\mu `$. In two dimensions, the ten states at the magnetizations $`M`$ and $`M+2`$ and energies $`E_0+4\eta `$ within the band of states to which the calculation was restricted, are connected by 12 transition variable pairs. Averaging over the 12 equations which determine $`\beta ^{}`$ and $`\mu `$ considerably increases the statistical relevance of the data. In three dimensions, the number of transition variable pairs is even larger.

Thus, the transition variable method directly provides the derivative $`\mu (E,M)`$ of the entropy which we need for the determination of its maximum.

From the zeros of $`\mu (E,M)`$ the spontaneous magnetization is determined with high precision as shown in Figure 1. The variation of the spontaneous magnetization per spin, $`m_{sp,N}=M_{sp,N}/N`$, as a function of the specific energy, $`\epsilon =E/N`$, is displayed in Figure 2 for two-dimensional Ising models on an infinite and on a small square lattice with $`32^2`$ spins. The curve for the infinite lattice has been drawn from the exact results for $`M(T)`$ and $`E(T)`$ . Both curves coincide at low energies. In our analysis, $`M_{sp}(E)`$ for energies lower than the critical value $`E_{c,\mathrm{}}`$ of the infinite system is fitted by

$$M_{sp}(E)=A\left|\epsilon ^{}\right|^{\beta _{\epsilon ,eff}(E)}$$
(9)
with an energy dependent effective exponent $`\beta _{\epsilon ,eff}(E)`$. Here, $`\epsilon ^{}=(E_{c,\mathrm{}}E)/(E_{c,\mathrm{}}E_g)`$ is the reduced energy. $`E_g`$ denotes the ground state energy. The effective exponent, which is given by $`\beta _{\epsilon ,eff}=\text{d}\mathrm{ln}M_{sp}/\text{d}\mathrm{ln}\left|\epsilon ^{}\right|`$, yields the critical exponent $`\beta _\epsilon `$ in the limit $`EE_{c,\mathrm{}}`$. Contrary to Ref. the critical expansion refers to $`E_{c,\mathrm{}}`$ and not to the critical energy of the finite system. The procedure has been performed for a series of system sizes $`N_1<N_2<N_3\mathrm{}`$ In the energy range where $`N_i`$ and $`N_{i+1}`$ yield the same exponent $`\beta _{\epsilon ,eff}(E)`$ this exponent is selected for the plot in Figure 3. This value of $`\beta _{\epsilon ,eff}(E)`$ is used up to the energy where the two exponents start to differ from each other, i.e. until finite size effects show up. Then we plot the common exponent for the system sizes $`N_{i+1}`$ and $`N_{i+2}`$ until these begin to disagree and so on. The largest systems for which we have performed these calculations were $`700\times 700`$ in two dimensions and $`80\times 80\times 80`$ in three (these are presumably the largest Ising systems for which microcanonical quantities have ever been determined). This approach has the advantage that the finite-size effects are monitored closely. Furthermore, one avoids to simulate unnecessarily large systems at low energy where small systems already yield data essentially free of finite-size effects.

The critical exponent $`\beta _\epsilon `$ is related to the commonly used exponent $`\beta `$ which describes the behavior of the order parameter with respect to the temperature by $`\beta _\epsilon =\beta /(1\alpha )`$ where $`\alpha `$ is the usual specific heat critical exponent . In two dimensions where $`\alpha =0`$ we expect $`\beta _\epsilon =1/8`$, in three dimensions with $`\beta =0.327`$ and $`\alpha =0.11`$ we expect $`\beta _\epsilon =0.367`$.

Also shown in Figure 3 is the corresponding effective exponent derived from the exact result $`M(E)`$ of the infinite two-dimensional Ising model. Obviously, the numerical data follow closely the exact result, thus demonstrating the reliabilty of our algorithm for computing accurate microcanonical quantities even for large systems.

Figure 3 exhibits two pecularities of the 2D system: (i) The strong $`N`$ dependence of $`E_{c,N}`$ makes it extremely difficult to obtain data points for $`\left|\epsilon ^{}\right|<0.2`$. (ii) The effective critical exponent $`\alpha _{eff}(t)`$ of the specific heat $`c(t)`$ of the infinite system is finite for finite $`t=(T_cT)/T_c`$, but precipitates to zero as $`t=0`$ is approached. This causes the dramatic behavior of $`\beta _{\epsilon ,eff}(\epsilon ^{})=\beta _{eff}/(1\alpha _{eff})`$ in the vicinity of $`\left|\epsilon ^{}\right|=0`$.

Remarkably, good values of the exponents are already found for the smallest systems considered here, namely $`L=50`$ in two and $`L=10`$ in three dimensions. This is especially true for the three-dimensional case where, over the whole energy range considered, the value of the effective exponent does not noticeably differ from the value of the critical exponent. Even for energies more than 30 % below the critical energy and for systems consisting of only 1000 spins, the obtained effective exponent is almost identical to the critical exponent. This is in marked contrast to a similar canonical analysis of the behaviour of the order parameter with respect to the temperature where large corrections to scaling lead to an effective temperature dependent exponent which differs by almost 50 % from the value of the critical exponent at temperatures some 30 % below $`T_c`$. Whereas it is obvious from Figure 3 that corrections are present when analysing the two-dimensional Ising model, these corrections are again small compared to the corrections observed for the temperature dependent effective exponent derived from the order parameter $`M(T)`$.

The origin of the different behaviour of $`\beta _{\epsilon ,eff}`$ as function of the reduced energy in two and three dimensions can be traced back to the logarithmic divergence of the specific heat in the two-dimensional Ising model. We expect a behaviour similar to that observed in the three-dimensional Ising model for all models with a non-vanishing specific heat critical exponent. Future investigations of various models will be needed to clarify this point.

## IV Summary

We have developed an algorithm which rapidly builds up transition variables for any reversible rule of jumping from one macrostate of a system, characterized by a set $`\underset{¯}{P}`$ of parameters, to another one. The most recent values of the transition variables govern the further path of the system through phase space in such a way that equally good statistics is obtained over the whole area which is covered by the calculation. As an example we have chosen the two- and three-dimensional Ising models with $`\underset{¯}{P}=(E,M)`$ together with single spin flips as the transition mechanism. From the transition observables we obtain highly accurate data for $`dS/dM`$ which allow a precise determination of the microcanonically defined spontaneous magnetization $`M_{sp}(E)`$. An analysis of $`M_{sp}(E)`$ for different system sizes yields excellent values for the order parameter critical exponent $`\beta _\epsilon `$.


format: text

text:
# Dense Antihydrogen: Its Production and Storage to Envision Antimatter Propulsion

## 1 Introduction

Within this century, the challenge to develop deep space missions will be undertaken. To complete a mission within a reasonable time frame, even to the nearest objects of interest, the Oort Cloud or the Alpha Centauri star system (4.3 light years away), the velocity of the spacecraft needs to be high, up to more than 10% of the speed of light. To achieve this one needs the highest energy-density fuel conceivable. This would be antimatter; a large amount of it and in a compact form.

Antimatter can produce three orders of magnitude more energy per gram than fission or fusion and ten orders of magnitude more energy than the chemical reactions currently used for propulsion. As a result, it is a prime candidate for use in future exploration beyond the solar system. It also is a candidate for future missions to the edge of the solar system, which now require on the order of 15 years after launch to reach Pluto.

An ideal antimatter propulsion system would convert all the energy from matter-antimatter annihilation into propulsion. Of course, it has long been realized that only a fraction of the total energy could be used in an antimatter propulsion system, but much thought has also gone into optimizing this fraction . Further, such a system would require what is by today’s standards a very large amount of stored antimatter. For a 500 kg probe, it would take more than 5 g of antimatter to get out to 200 Astronomical Units (past the Heliopause) in 1 year and more than 2.5 kg of antimatter to go 1 light year in 10 years.

Alternative propulsion schemes have also been proposed that use antimatter as a catalyst. They would require substantially less antimatter. (For example, they would use antimatter to drive fusion reactions  or fission reactions .) But for them the total weight of nuclear fuel for the same energy output is again roughly 3 orders of magnitude larger than for a pure annihilation energy source.

The amount of antimatter that can be stored per unit volume is critical to make the necessary equipment of reasonably low mass and to minimize the requirements on secondary technical support systems; i.e., cryogenic refrigerators. Unfortunately, the density of the most obvious candidate for a significant amount of stored antimatter, antiprotons, is limited by the Brillouin density ,

$$n_0=\frac{B^2}{8\pi mc^2}.$$
(1)
This is the limiting density for radial confinement in a magnetic field. For antiprotons stored in a 6 T magnetic field, typical for today’s technology, it would be around $`2\times 10^{12}`$ cm-3. Therefore, to achieve high storage densities, it will be necessary to make and store neutral antimatter. The only realistic neutral candidate is antihydrogen.

Further, a Bose-Einstein Condensate (BEC) of antihydrogen, even with present day BEC hydrogen densities, would be orders of magnitude more dense than the Brillouin storage density limit for charged antiprotons. This latter is itself orders of magnitude higher than the highest antiproton density so far achieved, $`10^6`$ cm-3 . To make a BEC of antihydrogen would be an individual step, where one could learn the techniques of controlling a relatively large amount of antihydrogen. The clear ultimate goal would be to make very dense antihydrogen, in the form of clusters or solids (perhaps stored diamagnetically).

Therefore, since a) controlled antiprotons have been created, b) cold antihydrogen has just been made, c) controlled antihydrogen is envisioned in the foreseeable future, and d) high-density BEC hydrogen has also been made (see the next paragraph), we find it appropriate to consider the possibility of producing antihydrogen BEC as a first, but critical, step towards the potential use of antimatter in a propulsion system.

Fundamental quantum theory predicts that there are ordinarily two kinds of particles; i) Fermions, which are particles with half-integer internal quantum angular momenta (like electrons and protons) and ii) Bosons, which are particles with integer internal quantum angular momenta (like photons or even numbers of Fermions such as the hydrogen atom). From the early work of Bose and Einstein, it was realized that many Bosons, if cool enough, could all fall into the lowest allowed quantum energy and behave as one single (even macroscopic) quantum object with remarkable properties. Two well-known examples are superconductivity, a sea of paired electrons, and superfluids, usually a liquid of 4He atoms (two protons, two neutrons, and two electrons). The recent discovery that atomic hydrogen can be condensed into a Bose-Einstein state, opened up new possibilities.

The purpose of this communication is to present an initial evaluation of possible routes towards efficient production and long-term storage of antimatter and also to identify the key issues that need to be addressed, specifically using BEC antihydrogen. We note that much of the equipment and the techniques that would be needed to study some of these concepts are already commonly available. Further, because of the symmetry between matter and antimatter, one should be able to perform many of the studies with matter. This is important since antiprotons are not (yet) readily available in large quantities.

## 2 Antihydrogen

Antihydrogen is the antimatter equivalent of hydrogen. It is composed of a positron in an atomic orbit around an antiproton. Antihydrogen is produced if a free antiproton ($`\overline{p}`$) picks up a positron ($`e^+`$) in an atomic orbit.

The antiprotons for antihydrogen can be made in high-energy collisions. They are currently being produced at the Fermilab and CERN particle accelerators for use in high-energy physics experiments. Antiprotons at CERN have been slowed and captured in Penning traps . The collection process is currently quite inefficient, but this efficiency could be improved by several orders of magnitude with straightforward changes . The positrons for antihydrogen are much more easily made using either small electron accelerators  or collected from certain radioactive sources with well-known gas moderation techniques .

Previously, a few atoms of antihydrogen have been produced at CERN  and Fermilab , but these have been produced at relativistic speeds, much too fast to capture. However, the ATHENA experiment at CERN recently produced a significant amount (of order 50,000) of cold antihydrogen atoms , a feat which caused international excitement . Later, the ATRAP collaboration also reported production of antihydrogen .

The antihydrogen program at CERN  is geared towards tests of fundamental physics principles such as spectroscopic studies of CPT invariance and the weak equivalence principle. While the aim of that research is to produce, what is for our purposes, a moderate number of anti-atoms, much can be learned from it about the different production and storage mechanisms proposed for our present purposes.

## 3 Production of Cold Antihydrogen

For any scheme proposed to produce slow antihydrogen, the velocity of the positron and antiproton relative to each other must be very small. This requirement is best described in terms of a temperature of the order of several Kelvin, describing the energy distribution in the center-of-mass frame. Such temperatures can be achieved by trapping the positrons and antiprotons in Penning traps. These traps confine charged particles in the radial direction with a strong solenoidal magnetic field (typically of the order of several Tesla), and in the longitudinal direction with static electric fields.

Once either positrons or electrons have been captured in a Penning trap, they cool very rapidly to equilibrium with the ambient temperature of the trap by emitting synchrotron radiation. In a 6 Tesla magnetic field the time constant for this process is approximately one hundred milliseconds.

Antiprotons have a much slower time constant for cooling through synchrotron radiation due to the heavier mass, but they can be cooled by loading electrons in the same trap. Coulomb collisions between the electrons and the antiprotons lead to thermalization of the two species, thus cooling the antiprotons. The heat input to the cold electrons due to this process is continuously removed by synchrotron radiation into the heat sink of the cryogenic environment of the experimental apparatus. Therefore, the velocities of all particles involved in this process can be reduced by cooling the physical structure of the trap to cryogenic temperatures.

This technique was first introduced by Rolston et al.  and has been successfully used by two experiments at CERN, PS196  and PS200 . The latter experiment currently holds the world record for the number of cold antiprotons confined in a trap of 1 million . The PS200 technology has now been fully incorporated into the ATHENA experiment at CERN .

## 4 Processes to Combine Antiprotons and Positrons into Antihydrogen

When an antiproton and a positron come together, they can only become an antihydrogen atom if there is a mechanism to remove energy from the system. This binding energy can be removed either by the radiation of a photon (radiative recombination) or it can be carried off by a third particle (three-body recombination). A number of different channels have been considered for such processes .

One of the strongest candidates for our purpose appears to be the three-body recombination (TBR) reaction $`\overline{p}+2e^+\overline{H}+e^+`$. The calculated rate for this reaction can be quite high :

$$\mathrm{\Gamma }=6\times 10^{13}\left(\frac{4.2}{T}\right)^{9/2}n_e^2[\mathrm{s}^1]$$
(2)
per antiproton, where $`n_e`$ is the positron density in number per cm3 and $`T`$ is the temperature. For a positron density of $`10^8`$ cm-3 at $`4.2^{}`$ K the half life for an antiproton in the positron plasma before it is converted into a antihydrogen atom is $`170\mu `$s. This essentially promises an instantaneous conversion of all antiprotons into antihydrogen atoms.

But most of these atoms are formed in highly excited Rydberg states, and it currently is not clear what will be the effect of the electric and magnetic fields in the trap environment on the net production rate. The recombination rate observed to date by the ATRAP collaboration is significantly lower than the predictions from equation (2). Most likely this is due to field ionization of the higher lying Rydberg states. To overcome this problem and to increase the net production rate it will be necessary to de-excite these atoms with an appropriate laser. Care must then be taken that photo-ionization of the Rydberg atoms does not limit the net production.

A second production mechanism of interest to us is the spontaneous radiative recombination (SRR) in which a positron is bound to an antiproton through emission of a photon. In this process 60% of the antihydrogen atoms are formed in states $`n=1`$ through $`10`$ and 20% of the atoms are formed directly in the ground state. The recombination coefficient $`\alpha ^{SRR}=\sigma ^{SRR}(v)v`$ for this process (at a positron temperature of 0.1 meV) is of the order of $`10^{10}`$ cm3/s. This recombination process is therefore substantially slower than the theoretical three-body reaction at the same low temperature. But it produces atoms stable against reionization by the strong electromagnetic fields in the environment. Additionally, the recombination rate can be enhanced by using a laser to stimulate a transition into a specific, stable, low-lying $`n`$ level.

Assuming the parameters proposed for the ATHENA experiment of $`10^7`$ antiprotons and $`10^8`$ positrons confined to a volume of 1 cm3, this reaction yields a production rate of 90,000 antihydrogen atoms per second. The entire sample of $`10^7`$ antiprotons would be converted to antihydrogen in less than 20 minutes. Considering the fact that the long range goal will be to store the produced antihydrogen for extended periods of time this rate is fully sufficient for initial experiments.

Several ideas to enhance the net production rate for both of the above mechanisms have been discussed in the literature .

## 5 Experimental Methods for Antihydrogen Production

Positrons and antiprotons have opposite charge, so they must be trapped in potential wells of opposite signs in Penning traps. Several ideas have been discussed in the literature on how to make the two oppositely charged clouds physically overlap in space, either in static configurations or by moving the electric wells through one another. The simplest rendition of such a scheme is the “nested trap” . In its most straight forward realization it consists of a sequence of five cylindrical electrodes that are biased in such a way that an elongated potential well is formed for antiprotons. This well which has a hill in its center that can then serve as a well for the oppositely charged positrons. (See Figure 1.)

To overcome the possible separation of the oppositely charged plasmas, the ATHENA collaboration uses an approach where the antiproton cloud is released from a side well into a positron plasma. The density and length of the positron plasma is chosen so that an individual low-energy antiproton entering the positron plasma along the axis will come to rest before leaving the plasma column again. Therefore, it will recombine with high probability in a single passage.

All of the above considerations have been successfully incorporated in the ATHENA experiment, which recently detected the creation of around 50,000 cold anti-hydrogen atoms at CERN . The creation was demonstrated by the detection of the annihilation products, preferentially on the walls of the trap. (See Figure 2.) This completes a major step on the road to the envisioned technology.

Once antihydrogen has been formed it will no longer be trapped in the Penning trap configuration. It will escape in the direction of the antiproton’s initial velocity and most likely will annihilate at the wall of the trap. This will make it difficult to diagnose the dynamics of the recombination reaction and the dependence of the recombination rate on the many different parameters relevant in these experiments. While recombination at rest will be an important part of an experiment aimed at filling some form of neutral storage scheme, we feel a method allowing easier diagnostics schemes is appropriate for the first stages of development.

Observe that it is the relative velocity between positrons and antiprotons which drives the reaction rate. The velocity of a positron at 4 Kelvin is about $`1.1\times 10^4`$ m/sec. An antiproton of the same velocity has an equivalent energy of about 0.6 eV (or 7300 K). This allows one to accelerate antiprotons into or through positron clouds without diminishing the conversion rate to antihydrogen.

Because of this, an appropriate proposal to make antihydrogen starts by keeping the antiprotons and positrons in separate potential wells in a Penning trap, as shown in Figure 3. The potential of the antiprotons can be raised with a small applied voltage so that when the antiprotons are released, they are accelerated along the axis of the trap. When they enter the positron plasma, some antiprotons will pick up positrons to form antihydrogen. The neutral antihydrogen will no longer be confined, so it will exit the trap along the axis, making a beam of antihydrogen.

For sufficiently large positron densities, the production efficiency should be very high. To produce a well defined beam, the initial velocity of the antiprotons and the density and length of the positron plasma must be carefully chosen.

## 6 Antihydrogen Storage

Long term storage of substantial amounts of antimatter must be developed to enable space missions relying on antimatter-based propulsion systems. Although it is clear that ultimately neutral antimatter must be used, up to now, no valid long-term storage concept for large quantities of antihydrogen has been developed. On the other hand, BEC confinement of neutral spin-polarized hydrogen atoms at densities up to $`5\times 10^{15}`$ cm-3 has been demonstrated. . We believe that developing the tools for making antihydrogen BEC could be an important step towards solving the antimatter storage problem.

While storage of high density samples of neutral alkali atoms has been demonstrated using magneto-optical and optical traps, the lack of appropriate laser systems has so far precluded hydrogen from such methods. Nevertheless, the BECs with the largest total number of particles have been formed using hydrogen. Until now neutral hydrogen has been stored either in cryogenic cells, by using wall collisions as the confinement mechanism, or in magnetic traps employing static magnetic gradient fields . Although total quantum reflection at containment walls is a possibility, without such a new technology, wall collisions are clearly unacceptable for the case of antihydrogen. One then needs to consider magnetic traps.

A static magnetic field trap is constructed with a minimum at its center. Then the only states that can be trapped are the so-called low-field seeking states. (The upper two states in the Hyperfine diagram for the ground state of atomic hydrogen.) These states will decay by the mechanism of dipole relaxation to the lower, un-trapped states with a rate constant of $`10^{15}cm^3s^1`$ . Another (smaller) loss mechanism is three-body recombination, which for hydrogen is ($`10^{19}cm^3s^1`$) . Apparently, these effects would strongly limit the density of achievable BECs.

But there are a number of other considerations that mean a controlled dense BEC is not ruled out in principle. Firstly, it is interesting to note that in a BEC, the coherence properties reduce the above dipole and three-body transition rates by factors of 2 and 6, respectively . Such a factor has been verified experimentally using Rubidium .

Next, Lovelace et al. showed  that a trap confining all hyperfine states of neutral (anti)hydrogen can be constructed using AC magnetic fields. (The situation is similar to strong focusing in particle accelerators.) Such a trap has not yet been experimentally realized because i) magnetostatic traps have been entirely sufficient for current research interests and ii) the technical challenges presented in generating AC magnetic fields with the appropriate amplitude and frequency are sizable. However, especially given the proper interest, traps of this type are certainly feasible.

But most importantly, the above transition rates are for current types of traps, with magnetic fields on the order of up to a few T. The rates are strongly magnetic-field dependent, and disappear for very large magnetic fields, of order 140 T . Fields on this order can now be made in a pulsed mode. Efforts to eventually create them statically have commercially reached 21 T .

At present the wasteful method of resonant evaporative cooling is used to achieve the temperatures and densities needed to form a hydrogen BEC. But the development of lasers for direct and efficient cooling of hydrogen atoms has now just started. Efficient laser cooling of hydrogen will revolutionize the methodology of forming, controlling, and studying hydrogen Bose-Einstein Condensates. All the cooling, as well as the optical storage, methods that have been so successful with Alkali BECs could then be used with antihydrogen. That is to say, once an antihydrogen atom has been made and stored, in principle the only additional requirement (beyond those needed for normal hydrogen) is the requirement that the residual pressure and temperature be low enough so that no residual normal atom be around to annihilate an antihydrogen.

Further, enormous advances have been achieved in the control of Alkali BECs leading to the demonstration of the coherent control of a BEC and the demonstration of the concept of atom-lasers . Based on these developments a concept of coherently amplified antimatter waves could be envisioned to increase antihydrogen storage capabilities by a significant factor. Such schemes are far beyond current technological capabilities, but the study of the underlying physical principles is a necessary beginning if one is to hope for a breakthrough discovery.

If the envisioned progress comes to fruition, laser cooling could then be used in an attempt to efficiently make an antihydrogen BEC.

## 7 Conclusion

An antihydrogen BEC would be an important step down a path that could eventually lead to even more dense anti-hydrogen molecules, liquids, solids, and cluster ions -. Indeed, since one might expect the next stage to be going from controlled ultra-cold (below 50 $`\mu `$K) BEC hydrogen atoms to controlled hydrogen molecules, it is heartening that there is evidence of a hydrogen-molecule superfluid with a critical temperature of 0.15 K . Since the triple point of hydrogen is at 13.8 K, a potential path to denser condensed antimatter becomes more interesting.

In our opinion a space-certified storage system for neutral antimatter can not be obtained from a linear extrapolation of heretofore existing technologies. Rather, it requires a scientific and/or technological breakthrough of the kind that Bose-Einstein Condensation might provide. While breakthroughs can never be predicted, they typically will not happen without the definition of a strong need and the challenge presented to the scientific community by a truly ambitious goal. Meanwhile many of the underlying issues can be addressed with both the modest supply of antimatter available at this time at accelerator centers world wide and with the limited means to store the particles. The technological and scientific knowledge gained in these tests will enable us to lay out a path into the future of antimatter-based propulsion systems.

Current antiproton production rates are low. While clever techniques can enhance these rates by several order of magnitude and quantities sufficient for advanced concepts can be produced given enough economic and political pressure onto the few available sources, a real breakthrough can only come through continued interest and research in this area. A good analogy is the comparison between a light bulb and a laser. In both cases light is produced, but in one system through thermal heating of a material and in the other through coherent processes. Antiprotons are currently produced by heating a metal target with a primary proton beam. This is a direct analogy to the light bulb — we are still awaiting the invention of a ‘laser-equivalent’ for the production of particles of antimatter.

Similarly, the efficient production and storage of antihydrogen (perhaps as a BEC, cluster, or solid), which are necessary steps towards large large-scale production as a fuel, can be studied with the quantities of antihydrogen soon to be available to us. The new modest production at CERN  should help us in identifying key problems and in developing technologies in preparation for future higher production rates of antimatter.

CERN is the only choice for the location to make the initial scientific advances. There are a number of steps needed to create and understand the desired form of antihydrogen. (i) Trap and cool the antihydrogen in a magnetic trap, developing tunable Lyman-alpha lasers along the way. After this first step one could (iiA) convert the antihydrogen into a BEC. Alternatively, one could (iiB) form and control molecular antihydrogen, and finally (iii) convert the antihydrogen into a superconducting cluster and/or solid of antihydrogen. Each of these advances, with strong funding, might be made in on the order of 5-10 years.

Future, even moderate, sources of antihydrogen (say at GSI or Fermilab) would also enable research on antimatter-based propulsion. One could also consider the use of antimatter to drive fusion reactions  or fission reactions  and also study crucial issues necessary to develop a propulsion system based entirely upon antimatter annihilation.

To meet the challenge of interstellar travel we have to overcome enormous scientific and technological barriers. Antimatter-matter annihilation is one of the prime candidates to achieve the high specific impulse i) desired for the challenging missions of exploring the Heliopause and visiting the Oort Cloud, and ii) needed if we plan to attempt a rendezvous with the nearest star systems. While no clear pathway to the necessary technologies exists, possible concepts can be identified based upon the need for the highest possible energy densities.

The experimental development of such systems as BECs and atom lasers in the normal matter world of laboratory-sized research equipment can help us to reach these most ambitious goals. To achieve them quickly it is necessary to set ourselves in motion now.

## Acknowledgements

We thank Michael Charlton, Lee Collins, Terry Goldman, TJ Johnson, Peter Milonni, Slava Turyshev, and especially Eddy Timmermans for their helpful comments and information. We all (MMN, MHH, and TJP) acknowledge the support of the US Department of Energy.

## References


text:
# Quenched Results for Light Quark Physics with Overlap Fermions Based on talks by L. Giusti, C. Hoelbling and C. Rebbi. Work supported in part under DOE grant DE-FG02-91ER40676.

## 1 Introduction

In the last few years a major breakthrough in the lattice regularization of Fermi fields was achieved through the closely related domain wall  and overlap  formulations. Neuberger found a lattice Dirac operator $`D`$  which avoids doubling and, most notably, satisfies the Ginsparg-Wilson (GW) relation . Thus the corresponding action in the massless limit preserves a lattice form of chiral symmetry at finite lattice spacing and volume . As a result the use of the overlap action entails many theoretical advantages (see ref.  and references therein for recent reviews); in particular it forbids mixing among operators of different chirality and therefore it can be very helpful (if not crucial) for computing weak amplitudes. No power-divergent subtractions are needed to calculate the matrix elements relevant for the $`\mathrm{\Delta }I=1/2`$ rule in $`K\pi \pi `$ decays .

The theoretical advantages of overlap fermions come at the cost of an increased computational burden in numerical simulations. The main question we wish to answer in the work presented here is whether the overlap formalism can be effectively used for large scale QCD calculations, with known algorithms and the current generation of computers. The calculation of light quark masses is an ideal test since it uses many of the ingredients needed for a generic phenomenological computation.

We have performed a fully non-perturbative calculation of $`(m_s+\widehat{m})`$ in the quenched approximation following the procedure proposed in . We also report results for the chiral condensate computed using the Gell-Mann–Oakes–Renner (GMOR) relation and for the pseudoscalar meson masses and decay constants.

## 2 Quark Masses and Chiral Condensate

The fermionic action in the overlap regularization reads

$$S=\overline{\psi }\left[\left(1\frac{1}{2\rho }aM\right)D+M\right]\psi $$
(1)
where $`M`$ is the matrix of bare masses $`(m_1,m_2,\mathrm{})`$ in flavor space. The Neuberger-Dirac operator is defined as

$`D`$ $`=`$ $`{\displaystyle \frac{\rho }{a}}\left(1+V\right)={\displaystyle \frac{\rho }{a}}\left(1+X{\displaystyle \frac{1}{\sqrt{X^{}X}}}\right)`$ (2)
$`X`$ $`=`$ $`D_W{\displaystyle \frac{1}{a}}\rho `$ (3)
where $`D_W`$ is the Wilson-Dirac operator (definitions and conventions used here are fully described in ref. ). The renormalized quark mass is defined as

$$\widehat{m}(\mu )=\underset{a0}{lim}Z_m(a\mu )m(a)$$
(4)
where $`Z_m(a\mu )`$ is a logarithmically divergent renormalization constant and $`m(a)`$ is the bare mass parameter which appears in the action of eq. (1) and in the corresponding axial and vector Ward identities. The non-singlet “local” bilinear quark operators we are interested in are defined as

$$O_\mathrm{\Gamma }(x)=\overline{\psi }_1(x)\mathrm{\Gamma }\left[(1\frac{a}{2\rho }D)\psi _2\right](x)$$
(5)
where $`O_\mathrm{\Gamma }\{V_\mu ,A_\mu ,S,P\}`$ correspond to $`\mathrm{\Gamma }\{\gamma _\mu ,\gamma _\mu \gamma _5,\text{11},\gamma _5\}`$. They are subject to multiplicative renormalization only, i.e. the corresponding renormalized operators are

$$\widehat{O}_\mathrm{\Gamma }(x,\mu )=\underset{a0}{lim}Z_\mathrm{\Gamma }(a\mu )O_\mathrm{\Gamma }(x,a)$$
(6)
where $`Z_\mathrm{\Gamma }(\mu )`$ are the appropriate renormalization constants. Since $`V_\mu `$, $`A_\mu `$ and $`S`$, $`P`$ belong to the same chiral multiplets $`Z_V=Z_A`$ and $`Z_S=Z_P`$. Flavor symmetry imposes $`Z_S=1/Z_m`$.

The bare chiral condensate is defined as

$$\chi (a)\underset{m0}{lim}\frac{1}{N_f}\overline{\psi }(0)[(1\frac{a}{2\rho }D)\psi ](0)$$
(7)
where $`m`$ in this case is a common mass given to the light quarks. For non-zero quark mass

$$\frac{1}{N_f}\overline{\psi }(0)[(1\frac{a}{2\rho }D)\psi ](0)=\chi (a)+\beta _\chi \frac{m(a)}{a^2}$$
(8)
since chiral symmetry forces the coefficient of the linear divergence to be zero. The condensate satisfies the integrated non-singlet chiral Ward identity

$`{\displaystyle \frac{1}{N_f}}\overline{\psi }(0)[(1{\displaystyle \frac{a}{2\rho }}D)\psi ](0)=`$ (9)
$`m{\displaystyle \underset{x}{}}P(x)\overline{P}(0)`$
Therefore by writing the correlation function $`P(x)\overline{P}(0)`$ as a time-ordered product and by inserting a complete set of states in standard fashion we can also write

$$\chi (a)=\underset{m0}{lim}\frac{m}{M_P^2}\left|0|P|P\right|^2$$
(10)
where $`M_P`$ is the mass of the pseudoscalar state $`|P`$. If we use

$$2m|0|P|\pi |=f_PM_P^2$$
(11)
where $`f_P`$ is the corresponding pseudoscalar decay constant, we arrive to the familiar GMOR relation

$$\chi (a)=\underset{m0}{lim}\frac{f_P^2M_P^2}{4m}$$
(12)
To preserve eq. (9), the renormalized chiral condensate is defined as

$$\frac{1}{N_f}\overline{\psi }\psi (\mu )=\underset{a0}{lim}Z_S(a\mu )\chi (a)$$
(13)
## 3 Non-perturbative Renormalization

The bilinear renormalization constants $`Z_\mathrm{\Gamma }(a\mu )`$ have been be computed at one loop in perturbation theory . Nevertheless, to avoid large uncertainties due to higher order terms, we prefer to compute them non-perturbatively.

The exact chiral symmetry implies a conserved axial current  which enters the corresponding Ward identities. Therefore the scheme and scale independent renormalization constant $`Z_A`$ can be computed from the axial Ward identity

$$Z_A\overline{}_\mu A_\mu (x)P(0)=2mP(x)P(0)+𝒪(a^2)$$
(14)
where $`\overline{}_\mu `$ is the symmetric lattice derivative. To compute the logarithmic divergent renormalization constant $`Z_S(\mu )`$ we implemented the RI/MOM non-perturbative renormalization technique proposed in ref. . The amputated off-shell Green’s functions are defined as

$$\mathrm{\Lambda }_\mathrm{\Gamma }(p)=𝒮^1(p)G_\mathrm{\Gamma }(p)𝒮^1(p)$$
(15)
where $`G_\mathrm{\Gamma }(p)`$ are the improved quark correlation functions and $`𝒮^1(p)`$ is the improved (external) quark propagator in the Fourier space. The projected amputated Green’s functions $`\mathrm{\Gamma }_\mathrm{\Gamma }(p)`$ are defined as

$$\mathrm{\Gamma }_\mathrm{\Gamma }(p)=\frac{1}{12}\text{Tr}\left[P_\mathrm{\Gamma }\mathrm{\Lambda }_\mathrm{\Gamma }(p)\right]$$
(16)
where the trace is over spin and color indices and $`P_\mathrm{\Gamma }`$ are the Dirac matrices which renders the tree-level values of $`\mathrm{\Gamma }_\mathrm{\Gamma }(p)`$ equal to unity (see  for more details). $`Z_S^{\text{RI}}(a\mu )`$ can be determined, up to $`𝒪(a^2)`$, by imposing the renormalization condition

$$Z_S^{\text{RI}}(a\mu )=\underset{m0}{lim}Z_A\frac{\mathrm{\Gamma }_A(p,m)}{\mathrm{\Gamma }_S(p,m)}|_{p^2=\mu ^2}$$
(17)
To compare the running of $`Z_S^{\text{RI}}(a\mu )`$ with the evolution predicted by the renormalization group equations, it is useful to define the renormalization group invariant (RGI) renormalization constant

$$Z_S^{\text{RGI}}(a)=\frac{Z_S^{\text{RI}}(a\mu )}{c_S^{\text{RI}}(\mu )}$$
(18)
where $`c_S^{\text{RI}}(\mu )`$ at the next-to-next-to-leading-order (NNLO) is given in . Up to higher order terms in continuum perturbation theory $`Z_S^{\text{RGI}}(a)`$ is independent of the renormalization scheme, of the external states and is gauge invariant. It is interesting to note the the non-degenerate scalar and pseudoscalar Green’s functions satisfy the WI

$`(m_1m_2)\mathrm{\Gamma }_S(p,m_1,m_2)=`$ (19)
$`m_1\mathrm{\Gamma }_P(p,m_1,m_1)m_2\mathrm{\Gamma }_P(p,m_2,m_2)`$
Equation (19) gives an exact relation among bare correlation functions of “local” scalar and pseudoscalar operators at non-zero quark masses, without any reference to the conserved axial or vector currents.

## 4 Numerical Results

The numerical results we present are based on a set of 54 quenched gauge configurations produced with the standard Wilson gauge action, with $`V=16^3\times 32`$ and $`\beta =6.0`$, which we retrieved from the “Gauge Connection” repository . For each configuration we fixed the Landau gauge by requiring a quality factor of $`\theta <10^6`$ and calculated the quark propagators for $`5`$ different bare quark masses $`m_i=\{0.040,0.055,0.070,0.085,0.0100\}`$, as described in appendix A. We computed then the two-point correlation functions

$`G_{SS}(t)`$ $`=`$ $`{\displaystyle \underset{x}{}}S(x,t)\overline{S}(0,0)`$ (20)
$`G_{PP}(t)`$ $`=`$ $`{\displaystyle \underset{x}{}}P(x,t)\overline{P}(0,0)`$ (21)
$`G_{AP}(t)`$ $`=`$ $`{\displaystyle \underset{x}{}}\overline{}_0A_0(x,t)\overline{P}(0,0)`$ (22)
We estimated the errors by a jackknife procedure, blocking the data in groups of three configurations, and we checked that blocking in groups of different size did not produce relevant changes in the error estimates.

### 4.1 Bare Quark Masses

In the quenched approximation the contributions from chiral zero modes is not suppressed by the fermionic determinant. In particular $`G_{PP}(t)`$ receives unsuppressed contributions proportional to $`1/m^2`$ and $`1/m`$, which should vanish in the infinite volume limit but can be quite sizeable at finite volume . These quenching artifacts cancel in the difference of pseudoscalar and scalar meson propagators

$$G_{PS}(t)=G_{PP}(t)G_{SS}(t)$$
(23)
The drawback is, of course, that the plateau in the effective mass and, correspondingly, the range that can be used for the $`\mathrm{cosh}`$ fit become shorter (See  for more details).

We fitted $`G_{PP}(t)`$ and $`G_{PS}(t)`$ to a single particle propagator in the time intervals $`t_1t_2=1216`$ and $`1016`$ respectively. The two correlation functions give values for the pseudoscalar masses compatible within the statistical errors. One does not notice any sign of the singular contributions from zero modes in the masses obtained from the pseudoscalar correlation function. These are expected to show up at some point, but one would probably need much higher statistical accuracy and lower values of $`m`$ to bring them into evidence. The results for the matrix elements, parameterised by the factors $`Z_{SP}`$ and $`Z_{PP}`$ for $`G_{SP}(t)`$ and $`G_{PP}(t)`$, respectively, show more significant differences. Within the large statistical errors, $`Z_{SP}`$ and $`Z_{PP}`$ are still compatible. However the central values are quite different and the fact that a curvature (see  for more details) shows up only in the results for $`Z_{PP}`$ points to the fact that what we are seeing is the effect of the unsuppressed zero modes, and not of chiral logarithms which would affect both sets of results (and would most likely become noticeable at much smaller values of $`am`$). On account of the above, we report in table 1 the results obtained from $`G_{SP}(t)`$ which we will use to derive our further results. It must also be said that most of the observables will be calculated directly at $`mm_s/2`$ (see below), and for these the difference between $`Z_{PP}`$ and $`Z_{SP}`$ is irrelevant within statistical errors.

We illustrate in fig. 1 the values for $`(aM_P)^2`$, obtained from $`G_{PS}(t)`$, as a function of the bare quark mass $`am`$. A linear behaviour

$$(aM_P)^2=𝒜_{M_P}+_{M_P}(am)$$
(24)
fits very well the data with

$$𝒜_{M_P}=0.0005(68)_{M_P}=1.43(7)$$
(25)
The vanishing of the intercept within statistical errors signals the absence of additive quark mass renormalization.

In order to fix the lattice spacing, we used the method of “lattice physical planes” , in which the ratio $`f_K/M_K`$ is fixed to its experimental value. This avoids recourse to a chiral extrapolation for observables except when it is really needed. For more details we refer the reader to  and only quote here the result

$$aM_K=0.216(8)af_K=0.0698(26)$$
(26)
To fix the lattice spacing we compared the value of $`af_K`$ with its experimental value and we obtained $`a_{f_K}^1=2.29(9)`$ We also computed the chiral condensate by using the GMOR relation. The quantity

$$a^3\chi _m=(am)\frac{Z_{PS}}{(aM_P)^2}$$
(27)
exhibits a very good linear behaviour in the bare quark mass and a linear fit leads to

$$a^3\chi =\mathrm{\hspace{0.17em}0.00117}(27)$$
(28)
The chiral condensate can also be computed directly from its definition in eq. (7). Using again a computational strategy similar to that used in eq. (23) in order to take care of the infrared divergent contributions (in $`m`$) from unsuppressed zero modes, we obtain

$$a^3\chi (a)=\mathrm{\hspace{0.17em}0.00117}(42)$$
(29)
which is in remarkable agreement with eq. (28).

### 4.2 Renormalization of the Axial Current

We implemented eq. (14) numerically by computing the ratio

$$R_\rho (t)=\frac{G_{AP}(t)}{G_{PP}(t)}$$
(30)
The results obtained for each mass are reported in the first plot of fig. 2. The flatness of the data highlights the improvement achieved with the overlap regularization. We fitted $`R_\rho (t)`$ to a constant $`\rho (m)`$ in the range $`5<t<27`$ and then performed a quadratic fit (second plot in fig. 2)

$$\rho (m)=A_\rho +B_\rho m+C_\rho m^2$$
(31)
obtaining

$`A_\rho `$ $`=`$ $`0.00002(7)B_\rho =1.286(3)`$ (32)
$`C_\rho `$ $`=`$ $`0.277(12)`$
Note that the intercept is compatible with zero: this should not come as a surprise since $`A_\mu (x)`$ has the correct behaviour under global non-singlet chiral transformations.

The axial current renormalization constant $`Z_A`$ is given by $`Z_A=\frac{2}{B_\rho }`$. We also fitted $`\rho (m)`$ linearly in the quark mass and we take the difference of the quadratic and linear results as a rough estimate of the systematic error. Our final result is $`Z_A=1.55(4)`$ . This value is larger than the one obtained in refs.  using standard lattice perturbation theory at one-loop.

### 4.3 Renormalization of the Scalar Density

To implement the RI/MOM non-perturbative renormalization technique we computed the projected, amputated Green’s functions of the quark bilinears defined in eq. (16). In the first plot of fig. 3 we show the ratio

$$R_{PS}=\frac{(m_1m_2)\mathrm{\Gamma }_S(p,m_1,m_2)}{m_1\mathrm{\Gamma }_P(p,m_1,m_1)m_2\mathrm{\Gamma }_P(p,m_2,m_2)}$$
(33)
for a given combination of masses. As expected from (19) its value is always compatible with one. The ratio of the vector and axial Green’s functions produces analogous results which we do not report here for lack of space. On the other hand the ratio $`\mathrm{\Gamma }_P/\mathrm{\Gamma }_S`$ with degenerate masses is expected to be sensitive to non-perturbative Goldstone pole contributions . A detailed analysis will be presented in a forthcoming paper.

To calculate $`Z_S`$ in the RI/MOM scheme, we computed the ratio

$$R_{AS}(m)=Z_A\frac{\mathrm{\Gamma }_A(p,m)}{\mathrm{\Gamma }_S(p,m)}$$
(34)
which we extrapolated to the chiral limit as

$$R_{AS}(m)=Z_S^{RI}+Am^2+Bm^2$$
(35)
to take into account the effect of unsuppressed zero modes. The results we thus obtained, shown in the second plot of fig. 3, lead to

$$Z_S^{\text{RI}}(2\text{GeV})=1.24(5)$$
(36)
where the error is mainly systematics and is due to the uncertainty in the value of the renormalization scale and to the chiral extrapolation (the latter estimated by using different forms of the extrapolating function). Using known results from NNLO continuum perturbation theory, we convert eq. (36) to the $`\overline{\text{MS}}`$ scheme, obtaining

$$Z_S^{\overline{\text{MS}}}(2\text{GeV})=1.41(6)$$
(37)
This result is in good agreement with the one in  which was determined by fixing the renormalized quark mass value. Invoking again NNLO continuum perturbation theory, we calculate $`Z_S^{RGI}`$ from eq. (18). We use $`N_f=0`$ and $`\mathrm{\Lambda }_{QCD}=0.238(19)`$ . The result, displayed in the second plot of fig. 3 shows the remarkable flatness of $`Z_S^{RGI}`$ for $`\mu >1.5`$ GeV which confirms the high level of improvement reached by the overlap regularization. Averaging over the region $`(1.53)`$ GeV, we find

$$Z_S^{RGI}=0.82(3)$$
(38)
## 5 Physical Results

In this section we will use our lattice results to infer the renormalized values of the sum of the strange and average up-down quark masses and of the chiral condensate. From eq. (26) we obtain

$$m_s(a)+\widehat{m}(a)=149(9)\text{MeV}$$
(39)
where the error is statistical only. Since our volume is fairly large, we expect our main sources of systematic errors to come from discretization effects and from the quenched approximation. For a rough estimate of the systematic error due to quenching approximation we can use the results in , from which one sees that, within the quenched approximation, the use of alternative observables to calculate the lattice spacing can produce differences of approximately $`10\%`$. Had we used $`r_0`$  to fix our lattice spacing we would have obtained a number $`7\%`$ higher than the one in eq. (39). Combining this fact with the results in , we infer that discretization uncertainties are below $`10\%`$. If we had used the experimental NNLO results from $`\alpha _s(M_Z)=0.118`$ in the matching of the renormalization constants in eq. (37), we would have obtained a value of the scalar renormalization constant $`10\%`$ higher than the one given above. The difference can be taken as an indication of the systematic error introduced by the quenched approximation. In order to be conservative, we will take $`15\%`$ as the estimate of our overall systematic error in the renormalized quark masses due to quenching and discretization effects. A more precise estimate of the systematic errors will need much more extensive simulations, which would go beyond the capability of our current computer resources and the exploratory nature of the our work. Combining the results in eqs. (36) and (39) we obtain

$$(m_s+\widehat{m})^{\text{RI}}(\text{2 GeV})=120\pm 7\pm 21\text{MeV}$$
(40)
which represent one of the main results of this work. By using the quark mass ratio $`m_s/\widehat{m}=24.4\pm 1.5`$ from chiral perturbation theory  the above translates to

$$m_s^{\overline{\text{MS}}}(\text{2 GeV})=102\pm 6\pm 18\text{MeV}$$
(41)
This result agrees very well with the current lattice world average .

Insofar as the value of the condensate is concerned, if we used the standard two-step approach, i.e. first measure the dimensionless condensate, see eqs. (28) or (29), and then multiply it by the cubic power of the lattice spacing, the result would be affected by a very large systematic error due to the uncertainty in the determination of the lattice spacing in quenched simulations. Instead, we will use an alternative method . We write the GMOR relation (12) for the renormalized condensate as follows

$$\chi (a)=\frac{1}{4}f_\chi ^2_{M_P}a^1$$
(42)
where $`_{M_P}`$ is defined in eq. (25) and $`f_\chi =0.1282`$ GeV is the “experimental” value, in physical units, of the pseudoscalar decay constant extrapolated to the chiral limit. While computing the condensate from the above formula relies on an additional element of experimental information, it has several advantages. The most important is that, by expressing the condensate in terms of $`f_\chi `$, we are left with only one power of the UV cutoff $`a^1`$. With this method we obtain

$`\overline{\psi }\psi ^{\text{RI}}(\text{2 GeV})=`$ (43)
$`\mathrm{\hspace{0.17em}0.0167}\pm 0.0010\pm 0.0029\text{GeV}^3`$
where the estimate of the systematic error has been made using the same criteria we used to compute the error for the quark masses. This is our best value for the chiral condensate. It is interesting to note that, if we had used the standard technique, starting from eq. (28), we would have obtained a result with a central value very close to the value in eq. (43), but with a much higher error. Finally, from eq. (43) and NNLO matching, we get

$`\overline{\psi }\psi ^{\overline{\text{MS}}}(\text{2 GeV})=\left(267\pm 5\pm 15\text{MeV}\right)^3`$ (44)
This result is in very good agreement with the result obtained by the authors of refs. , while it is smaller than the result in , even if still compatible within errors. Our result is also compatible within errors with the number obtained few years ago in  with Wilson-type fermions. We expect, though, the systematics due to the discretization effects to be smaller ($`𝒪(a^2)`$) in the result reported in eq. (44) than the error ($`𝒪(a)`$) which affects the determination in .

## 6 Conclusions

At this conference we presented results for the pseudoscalar meson masses and decay constants together with the first fully non-perturbative computation of $`(m_s+\widehat{m})`$ with overlap fermions in the quenched approximation. We also computed the chiral condensate $`\overline{\psi }\psi `$ from the GMOR relation and directly. To avoid uncertainties due to lattice perturbation theory, we computed the multiplicative renormalization constant $`Z_S(\mu a)`$ non-perturbatively in the RI/MOM scheme. While the systematics errors due to quenching are common to previous calculations, the other systematic errors (mostly discretization effects) are different than in other lattice regularizations and likely to be smaller, because of chiral symmetry. Our results have indeed produced a remarkable verification of “good chiral behavior” of the overlap fermions both in the axial Ward identity and for the pseudoscalar masses.

The calculation of light quark masses uses many of the ingredients needed for a lattice calculation of weak matrix elements, although the latter is computationally more demanding. From this point of view, the very good agreement between our results for the quark masses and the current lattice world average bodes well for the use of the overlap formalism also in matrix element calculations. Our investigation has been mostly of exploratory nature. One would need to extend it to larger volumes and better statistics. Nevertheless, we believe that it demonstrated that the overlap formalism can be used effectively, with known algorithms and the present generation of computers, for large scale QCD calculations, at least in the quenched approximation. Thus we would conclude that it represents a very promising non-perturbative regularization for solving long standing problems, such as the proof of the $`\mathrm{\Delta }I=1/2`$ rule and the calculation of $`ϵ^{}/ϵ`$.

## Appendix A Implementation of the Overlap Operator

The complicated form of Neuberger operator renders its numerical implementation more involved and expensive than those of the most common regularizations. The crucial point is to implement the sign function

$$ϵ(H)=\text{sign}(H)$$
(45)
over the whole spectrum of eigenvalues of $`H=\gamma _5X`$. Exact diagonalization is feasible  for two dimensional models, but becomes prohibitive for larger systems such as QCD. Many algorithms have been proposed in the literature . We opted for the optimal rational function approximation suggested in refs. . Starting from the observation  that the sign function can be written as

$`\epsilon (x)`$ $`=`$ $`\underset{N\mathrm{}}{lim}\epsilon _N(x)`$
$`\epsilon _N(x)`$ $`=`$ $`x\left(c_0+{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \frac{c_k}{x^2+q_k}}\right)`$ (46)
a good approximation can be obtained at finite $`N`$ by fixing the coefficients $`c_k`$ and $`q_k`$ with the Remez algorithm .

The numerical procedure we have followed to obtain a quark propagator from a fixed source vector $`|\eta `$, i.e. solving the equation

$$\left[(1\frac{1}{2\rho }m_i)D+m_i\right]|\chi _i=|\eta $$
(47)
is the following:

* We have extracted the $`15`$ smallest eigenvalues of $`H^2`$ and the corresponding eigenvectors minimizing the Ritz functional , projected them out and treated the action of $`ϵ(H)`$ on this subspace exactly. The minimum precision required for each eigenvalue $`\lambda _i^2`$ is

$$r_{ritz}\sqrt{\frac{|H^2|\lambda _i\lambda _i^2|\lambda _i|^2}{\lambda _i|\lambda _i}}<10^6$$
(48)
where $`|\lambda _i`$ is the corresponding eigenvector. We checked on all the configurations that the largest eigenvalue extracted is always $`|\lambda _{max}|>0.15`$.
* We have rescaled the reduced matrix $`\stackrel{~}{H}=H/r`$ by a factor $`r=3.7`$ and we have approximated $`ϵ(\stackrel{~}{H})`$ as in eq. (46) with $`N=14`$ and $`c_k,q_k`$ fixed with the Remez algorithm  which guarantees a precision

$$|ϵ(\stackrel{~}{H})ϵ_N(\stackrel{~}{H})|3.6\times 10^6$$
(49)
in the range $`\stackrel{~}{H}[0.040,1.8]`$. In order to invert all the equations

$$\left(\stackrel{~}{H}^2+q_k\right)|\psi _k=|sN=1,\mathrm{}14$$
(50)
in one stroke we use an (inner) multi conjugate gradient (MCG) solver . Note, that during the inner MCG it is necessary to store only $`N+2`$ large vectors since we are only interested in

$$|\psi =\underset{k=1}{\overset{N}{}}c_k|\psi _k$$
(51)
The convergence is governed by the inversion corresponding to the smallest $`q_k`$ for which we require a residual

$$r_{in}\sqrt{\frac{|(\stackrel{~}{H}^2+q_{14})|\psi _{14}|s|^2}{s|s}}<10^6$$
(52)
Notice that having extracted the lowest eigenvectors of $`H^2`$ improves the approximation used for $`\epsilon (H)`$ and reduces the condition number of $`H^2`$, which results in a speed up of the inversion in eq. (50).
* From eq. (47) and using the GWR we can write

$$|\chi _i=\frac{\rho \left((1\frac{m_i}{2\rho })D^{}+m_i\right)}{\left(\rho ^2m_i^2/4\right)\left(D+D^{}\right)+\rho m_i^2}|\eta $$
(53)
and therefore we are left to invert

$$\left(D+D^{}+\frac{\rho m_i^2}{\rho ^2m_i^2/4}\right)|\stackrel{~}{\chi }_i=|\eta $$
(54)
$`D+D^{}`$ is an Hermitian matrix and therefore a MCG can be applied also in this case. Moreover $`[D+D^+,\gamma _5]=0`$ implies that we can always use sources $`|\eta `$ and solutions $`|\chi _i`$ restricted to one chiral sector, saving one application of $`\epsilon (H)`$ per iteration . The convergence is governed by the inversion corresponding to the smallest quark mass, for which we require a residual of $`r_{out}<10^5`$

We did not attempt to use any preconditioning for the inner and the outer MCG. At the $`i`$-th step of the outer MCG, the difference $`|d_i`$ of two consecutive source vectors $`|s_i`$ and $`|s_{i1}`$ on which $`\epsilon _N(H)`$ has to be applied gets smaller and smaller when $`r_{out}`$ decreases. Therefore at the $`i+1`$ step we have applied $`\epsilon _N(H)`$ directly on $`|d_i`$ with an increased $`r_{in}`$ (every $`40`$ outer MCG step we required always the full precision to avoid rounding accumulation). After convergence, we always checked that the true residual associated with the eq. (47) for each mass is $`<10^5`$. In our simulations the average numbers of inner and outer iterations is $`𝒪(150)`$ and $`𝒪(250)`$ respectively.


text:
# 1 Introduction

## 1 Introduction

The Soft Gamma Repeaters are young neutron stars detected through their brilliant flashes of X-rays and gamma-rays. Strong physical arguments point to the presence of $`10^{14}10^{15}`$ G magnetic fields in these sources. (A detailed analysis of 27 August 1998 flare data is given in Feroci et al. 2001, and of the flare physics in Thompson & Duncan 2001.) The Anomalous X-ray Pulsars have gradually emerged as a class of neutron stars with moderate luminosities ($`L_X3\times 10^{34}10^{36}`$ erg s-1) and a narrow range of spin periods ($`P612`$ s). The AXPs are persistently spinning down, some are observed to vary in brightness, and about half are convincingly associated with young supernova remants. Remarkably, the SGRs in their quiescent states have similar luminosities and spins to the AXPs, and are spinning down even more rapidly. This overlap between the two classes of sources within a three-dimensional parameter space (in spite of the very different detection methods) motivated the suggestion that both are isolated neutron stars powered by the same energy source: the decay of an ultrastrong magnetic field (Thompson & Duncan 1993, 1996). Recent data suggest two objects as possible ‘missing links’ between the two classes: SGR 0525-66, which has been dormant as a burst source since 1983 and has a soft AXP-like power-law spectrum (Kulkarni et al. 2000); and the anomalous pulsar 1E 1048.1-5937, which has a harder persistent X-ray spectrum than SGR 0525-66, and a noisier spindown torque than the other AXPs (Kaspi et al. 2000). Nonetheless, the observed properties of the AXPs are less constraining than those of the SGRs, and have inspired a number of models over the last decade.

## 2 Models of the Anomalous X-ray Pulsars

The AXPs were first defined as a group111The name ‘Anomalous X-ray Pulsar’ was coined by Duncan & Thompson (1995) in a conference presentation. by Mereghetti & Stella (1995). Theoretical models for the AXPs must explain several peculiar properties of these sources: i) the absence of a modulation of the X-ray pulse frequency by binary motion, which is consistent with a low mass helium star, or a very low mass degenerate dwarf companion (Mereghetti, Israel, & Stella 1998); ii) low optical fluxes, $`L_{\mathrm{opt}}/L_X10^3`$ (Hulleman et al. 2000a,b); iii) persistent spindown, $`\dot{P}>0`$, much larger than the $`\dot{P}`$ of radio pulsars, and young characteristic ages $`P/\dot{P}=4\times 10^34\times 10^5`$ yrs; iv) a wide range of torque noise, with some AXPs spinning down almost as quietly as radio pulsars (Kaspi et al. 1999, 2001); iv) spindown luminosities $`I\mathrm{\Omega }\dot{\mathrm{\Omega }}10^{32}10^{33}`$ erg s-1, some 2-3 orders of magnitude smaller than the X-ray luminosities $`L_X3\times 10^{34}10^{36}`$ erg s-1; and v) the localisation of at least 3 AXPs close to the geometrical centers of young supernova remnants (e.g. Gaensler et al. 2001). See the paper by Israel, Mereghetti, & Stella (these proceedings) for a more complete review.

Here is an inventory of the models for the AXPs, listed in the chronological order in which they were invented. I will indulge in one editorial comment at the beginning: In constructing models, there has sometimes been a tendency to equate the ‘symptoms’ with the underlying ‘disease’. Before the splendid phase-connected spin measurements by Kaspi et al. (of the AXPs) and Woods et al. (of the the SGRs) the only neutron stars showing large torque noise were accreting X-ray binaries and a few young glitching pulsars (such as Vela). The most conservative approach would be to fit the AXP/SGR measurements into one of these classes, and this indeed has been the strongest motivation for the accretion-based torque models. There is, however, a fairly direct argument that we are observing a new mechanism for torque noise in neutron stars. The torque acting on an isolated neutron star is a direct measure of the current flowing through its outer magnetosphere; but the magnetic fields of the SGR sources are time variable (giving rise to the bright X-ray bursts), and so one can also expect the current to be time-variable. The only leap one must make here, with regard to the AXP sources, is to postulate that the magnetic field (and the crust which anchors it) are capable of more gradual deformations, in addition to the sudden yields and fractures which appear to trigger SGR flares. Nonetheless, the patterns of torque noise in the SGRs and AXPs are turning out to be sufficiently similar, that any torque model should be able to accomodate both classes of sources.

i) Massive white dwarfs, formed in white dwarf mergers, now spinning down by magnetic dipole torques. Paczyński (1990) suggested this model for 1E 2259+586, which was recognized early on as a peculiar object (Fahlman & Gregory 1981). The spindown luminosity of a white dwarf is much larger, at fixed $`P`$ and $`\dot{P}`$, than that of a neutron star, by the ratio of moments of inertia $`I_{\mathrm{wd}}/I_{\mathrm{ns}}0.3(M_{\mathrm{wd}}/M_{\mathrm{ns}})(R_{\mathrm{wd}}/R_{\mathrm{ns}})^2>3\times 10^4`$. The required polar magnetic field is

$$B_{\mathrm{pole}}R_{}^2=\frac{9\times 10^{26}}{\mathrm{sin}\alpha }\left(\frac{P\dot{P}}{10^{10}\mathrm{s}}\right)^{1/2}\left(\frac{M_{}}{M_{}}\right)^{1/2}\left(\frac{I_{}}{M_{}R_{}^2}\right)^{1/2}\mathrm{G}\mathrm{cm}^2,$$
(1)
from the standard dipole formula with inclination angle $`\alpha `$. King et al. (2001) have noted that a binary merger is a promising formation scenario for the white dwarf RE J0317-853 ($`B48\times 10^8`$ G), but its $`725`$-s spin is very slow for the AXPs (albiet fast for a white dwarf). The questions which hang over the white dwarf interpretation of the AXPs are fairly obvious: why is the radiative area deduced from the thermal component of AXP emission comparable to that of a neutron star (Perna et al. 2001)? Why are half the AXPs situated close to the centers of supernova remnants, whose kinetic energy greatly exceeds the rotational energy of the white dwarf? Why are the white dwarfs visible only while young, and why the narrow range of spin periods?

ii) Magnetars: neutron stars with decaying $`10^{14}10^{15}`$ G magnetic fields. Thompson & Duncan (1993, 1996) drew a connection between 1E 2259+586 and the SGR 0525-66, based on their similar spin periods, persistent $`L_X`$, and residence in $`10^4`$-yr old SNR. The rapid spindown of the AXPs has a simple interpretation as magnetic dipole breaking, modified by variable magnetospheric currents flowing along the decaying magnetic field (Thompson, Lyutikov, & Kulkarni 2001). A variable X-ray flux also points to active field decay in some sources. The most interesting argument against accretion as the energy source of AXP emission is indirect and comes from the SGR sources: in particular, SGR 1900+14 was detected as a persistent source within $`10^3`$ s of the August 27 giant flare, and yet the radiative momentum of the flare would excavate any accretion disk and suppress accretion for a much longer period (Thompson et al. 2000).

The limited duration of AXP/SGR activity is consistent with the freeze-out of ambipolar diffusion of the magnetic field through the core of a neutron star after surface photon cooling dominates (Thompson & Duncan 1996; Heyl & Kulkarni 1998). The characteristic timescale is $`10^5`$ yrs for an iron surface composition, and $`10^4`$ yrs for H/He. However, the observed range of spin periods is surprisingly narrow, given the wide range of characteristic ages exhibited by magnetar candidates: from $`P/\dot{P}=13,000`$ yrs for the spinning down SGRs 1806-20 and 1900+14 (Kouveliotou et al. 1998; 1999) to $`P/\dot{P}=4\times 10^5`$ yrs for 1E 2259+586. This points, in a model-independent manner, to a decay in the torque acting on 1E 2259+586, from its historic average (Thompson et al. 2000). The measured spins are remarkably close to the upward extrapolation of the radio pulsar death line in the $`\dot{P}P`$ plane. If the magnetar model were to explain them, it would probably be on the basis of a connection between a pair-loaded current, and the stability of a non-potential magnetic field outside the star.

iii) Low Mass X-ray Binaries. In defining the AXPs, Mereghetti & Stella (1985) suggested that they are accreting neutron stars with very low mass companions. They were motivated by the similarity between the spin periods of the AXPs and the equilibrium spin period $`10(B/3\times 10^{11}\mathrm{G})^{6/7}(L_X/10^{35}\mathrm{erg}\mathrm{s}^1)^{3/7}`$ s of an accreting neutron star whose magnetic field is typical of (in fact, slightly weaker than) a young radio pulsar. In addition, the known LMXB 4U 1626-67 has a spin and X-ray luminosity ($`P=7.7`$ s, $`L_X10^{36}`$ erg s-1) in the AXP range. The X-ray spectrum of 4U 1626-67 is, however, much harder than that of any AXP, and there is an optical counterpart (Chakrabarty 1998). This system also has a feature which is natural in any LMXB model for the AXPs, but is not observed in those sources as presently defined: 4U 1626-67 has exhibited intervals of rapid spin-up as well as spin-down, which is expected as the source is old. The long evolutionary timescale leads more generally to a prediction of too many AXP sources, given the probable association of 3 AXPs with $`10^4`$-yr old supernova remnants: a binary with a companion of mass $`M_20.1M_{}`$, undergoing Roche-lobe overflow with a luminosity $`L_X=GM_{\mathrm{NS}}\dot{M}/R_{\mathrm{NS}}=10^{35}`$ erg s-1, would evolve on the much longer timescale $`M_2/\dot{M}1\times 10^{10}(M_2/0.1M_{})(L_X/10^{35}\mathrm{erg}\mathrm{s}^1)^1`$ yrs. Overall, it is not surprising that at least one of the large number of accreting X-ray binaries should happen to overlap the AXP parameter space.

iv) Neutron Stars with Fossil Disks. Some neutron stars may capture a quasi-Keplerian disk after their formation, a scenario which has received renewed interested as a model for the AXPs (Van Paradijs, Taam, & Van den Heuvel 1995; Chatterjee, Hernquist, & Narayan 2000). As was realized after the discovery of the first pulsar planetary system, several evolutionary pathways could result in the formation of such a disk. The accretion rate will generically drop off with time, as the orbiting material spreads outward. If the spin of the neutron star is able to track the equilibrium spin period $`P_{\mathrm{eq}}\dot{M}^{3/7}`$ (where the corotation radius is close to the Alfvén radius) then its spin will decrease with time. Li (1999) has noted that in some sources (e.g. 1E 1048.1-5937) the time to return to the equilibrium spin is much longer than the measured characteristic age $`P/2\dot{P}`$, unless the X-ray emission is significantly beamed and $`\dot{M}`$ is underestimated. This is not a concern in the case of an old accreting binary with fluctuating $`\dot{P}`$, but it is if the source is only $`10^4`$ years old.

Several questions hang over this model: i) Why are the AXPs so similar to the SGRs in their quiesecent states? It has been suggested that the torque is modified by accretion effects in the SGR sources (Marsden et al. 2001), but no cogent connection between accretion and bursting activity has been given. ii) An upper age limit of $`10^410^5`$ yrs could result from the dissipation of the disk (e.g. Chatterjee & Hernquist 2000), but a narrow range of $`P`$ would seem to require a narrow range of initial disk masses. iii) Optical/IR emission from the portion of the disk exposed to the central X-ray source. The possible detection of an optical counterpart to the AXP 4U 0142+61 with (de-reddened) flux $`10^3`$ of the X-ray flux (Hulleman et al. 2000b) conflicts with the expectation of $`L_{\mathrm{opt}}/L_X10^2`$ (Perna, Hernquist, & Narayan 2000), and the spectrum is much redder than is typical of a LMXB. iv) A basic theoretical question is whether an isolated neutron star, surrounded by a centrifugally supported envelope, is able to lose angular momentum through a propeller and approach corotation with a weak accretion flow. Whether this is possible depends on the angular momentum carried away by each gram of ejected material, the initial spin, and the disposition of the surrounding material (in a disk vs. an extended envelope). Calculations which answer in the affirmative have assumed a disk geometry and a fast initial spin, but have allowed the ejected material to be spun up to corotation with the star – the most optimistic prescription for the propeller (Chatterjee et al. 2000). The critical dipole field, above which the rotation is slowed down to allow accretion, rises from $`10^{13}`$ G to more than $`10^{15}`$ G if one makes a different (weaker) prescription for the propeller: that the angular momentum per gram is that of a Keplerian orbit at the magnetospheric boundary.

## 3 Magnetars: Recent Developments

The magnetar model for the Soft Gamma Repeater sources and the Anomalous X-ray Pulsars has been reviewed recently (Thompson 2000), and so I will list a few recent theoretical highlights here. A case can be made that we now have a better understanding of the radiative mechanism and the nature of the ‘engine’ in SGR bursts, than has yet been obtained for the more distant cosmological GRBs.

i) Direct Evidence for Trapped Fireballs in SGR Giant Flares. After the initial $`0.5`$-s fireball, followed by a $`40`$-s phase of irregular pulsations, the X-ray flux in the 27 August giant flare maintained a smooth decline until an abrupt termination at $`380`$ s. After averaging over the large-amplitude $`5.16`$-s pulsations, the light curve is well fit by the contracting surface of a trapped fireball, $`L_X(t)L_X(0)(1t/t_{\mathrm{evap}})^3`$ (Feroci et al. 2001). Because the shape of the lightcurve is changed by neutrino cooling (which is strongly temperature-dependent), the diameter of the confined $`e^\pm `$ plasma cannot be much less than 10 km. This leads to a strict lower bound $`m/R_{\mathrm{NS}}^33\times 10^{13}`$ G on the dipole magnetic moment $`m`$ of the confining field (Thompson & Duncan 2001). As the X-ray flux declines, the thermal component of the X-ray spectrum maintains a constant temperature $`11`$ keV (black body value: Mazets et al. 1999; Feroci et al. 2000), identical to the temperature at which photon splitting freezes out (Thompson & Duncan 1995).

Independent evidence for a trapped fireball comes from a less energetic burst emitted by SGR 1900+14 two days after the 27 August giant flare. The main burst was followed by a faint pulsating tail lasting at least $`10^3`$ s, during which the excess flux decreases as $`t^{0.6}`$ (Ibrahim et al. 2001). The high temperature of this tail implies a radiative area only $`1`$ percent that of a neutron star. This area is consistent with a trapped fireball of a small volume, but an energy density similar to that attained in the much longer pulsating tail of the giant flare. The energy released in this faint tail can be accounted for by heating of the outer crust through $`e^{}`$-captures, while it was compressed under the fireball pressure (Ibrahim et al. 2001).

ii) Thermal Stability of Magnetized Pair Plasmas, and Short SGR Bursts. The trapped fireball model requires that energy is released fast enough into the magnetosphere to create a plasma in local thermodynamic equilibrium. A steady balance between heating and radiative cooling from the (optically thick) plasma could be maintained if the rate of injection of energy were less than $`10^{42}(R/10\mathrm{km})`$ erg s-1 within a volume $`R^3`$ (Thompson & Duncan 2001). We know that LTE can be achieved in the giant flares, because their peak luminosities exceed $`10^{44}`$ erg s-1 during the initial $`0.5`$-s transient. Similarly, one infers that the injection was fast enough to reach LTE in the main, bright component of the 29 August burst. Nonetheless, the question remains open as to whether an LTE fireball actually forms in the much more common short SGR bursts, which have a characteristic $`0.1`$-s duration much shorter than the $`3.5`$-s duration of the 29 August burst (Göğüş et al. 2001). In a steady state, the heated electrons cool primarily by Compton upscattering O-mode photons (which maintain a Wien distribution and have a large scattering cross-section $`\sigma _T`$ in an ultra-strong magnetic field). The net result is an anti-correlation between hardness and luminosity (Thompson & Duncan 2001) which is qualitatively similar to that observed in the short bursts (Göğüş et al. 2001).

iii) Passively Cooling, Ultramagnetized Neutron Stars. The transparency of the outer envelope of a neutron star is increased slightly in the presence of a strong magnetic field, which led to the suggestion that some AXPs may be passively cooling neutron stars without active magnetic field transport (Heyl & Hernquist 1997). This model is simple enough for detailed calculations of the angular and spectral distributions of the cooling X-ray flux to be compared with AXP data. Gravitational bending of the photon trajectories tends to reduce pulsed fractions (Psaltis, Özel, & DeDeo 2000). In addition, the spectrum will differ from a pure blackbody, as the result of the strong asymmetry between the opacities of the two X-ray polarization modes, which might explain some of the soft, high energy power-law tails obtained in AXP spectral fits (Özel 2001). It is worth keeping in mind that two effects will bias the number of observed sources toward those with active magnetic field decay: external currents are more efficient at generating X-rays than is deep heating; and the lifetime of a magnetar as a thermal X-ray source is lengthened, by an order of magnitude, by the decaying field (Thompson & Duncan 1995; Heyl & Kulkarni 1998).

iv) Magnetar Electrodynamics. The persistent emission of the active SGRs has a hard, power-law spectrum and correlates directly with X-ray burst activity (Woods et al. 2001). A strong magnetic field will heat the interior of a neutron star, and also drive electrical currents through its exterior (Thompson & Duncan 1996). The lowest energy deformation of the crust, extending over a scale of kilometers, involves a predominantly rotational motion, which will twist up the external magnetic field. Detailed force-free models of twisted, current-carrying magnetospheres have been constructed by Thompson, Lyutikov & Kulkarni (2001). They have the remarkable property that the optical depth to resonant cyclotron scattering depends directly on the twist angle, but not on the mass, charge, or radius of the scattering particle. As a result, the smooth pulse profile emerging in SGR 1900+14 after the 27 August giant flare (Woods et al. 2001) has a simple explanation through the formation of an extended scattering screen at the electron cyclotron resonance, which for electrons sits at $`R/R_{\mathrm{NS}}10(B_{\mathrm{pole}}/10^{14}\mathrm{G})^{1/3}(\mathrm{}\omega /\mathrm{keV})^{1/3}`$ in magnetar-strength fields. In addition, multiple resonant scattering will generate an extended high-energy tail in the persistent X-ray spectrum, and allow large pulse fractions in AXP sources. The field decreases more slowly with radius than a pure dipole, but only slightly, so that the polar surface fields inferred from the measured $`\dot{P}`$ and $`P`$ of the SGR sources are reduced by a factor $`310`$ compared with respect to the classical magnetic dipole model.

## 4 Observational Tests and their Theoretical Implications

In the next few years, the following observational efforts can be expected to throw further light on the true nature of the SGR/AXP sources.

i) Optical counterparts to the SGRs and AXPs. Besides the possible detection $`L_{\mathrm{opt}}/L_X10^3`$ for 4U 0142+61 (with a peculiar spectrum: Hulleman et al. 2000), there is an interesting upper limit for 1E 2259+586 (Hulleman et al. 2000a) and SGR 0525-66 (Kaplan et al. 2001). High extinctions may limit further progress on the SGR sources, but the prospects are good for improving optical identifications of the AXPs by making use of improved X-ray localizations. Strong constraints on accretion will be obtained – either through flux limits; or from the spectral distribution of the measured optical/IR emission, combined with the cross-correlation between the X-ray and optical fluxes.

ii) Connection between the SGR/AXPs and nearby cooling neutron stars. The SGRs and AXPs fade as X-ray sources after $`10^410^5`$ yrs, but little is known about the further evolution of their X-ray flux and spectrum. The nearby soft X-ray pulsars RX J0420.0-5022, RX J0720.4-3125, and RBS1223 have similar spin periods (22.7, 8.37, and 5.2 s) to the SGRs/AXPs, but much lower X-ray luminosities. It is difficult to make unambiguous predictions about the field decay occuring in a magnetar after it begins photon cooling, because the heating quickly becomes dominated by the crustal magnetic field. Measurements of spindown in these sources will provide crucial input to the physical modelling.

iii) Ion cyclotron features: emission vs. absorption? Detailed calculations of cooling neutron star atmospheres in very strong magnetic fields ($`B=10^{14}10^{15}`$ G) show that broad absorption features are possible at the fundamental ion cyclotron resonance (Zane et al. 2001; Ho & Lai 2001). Nonetheless, the persistent emission of the SGRs and some AXPs has a strong non-thermal component which, if interpreted in terms of stationary magnetospheric currents, implies a significant amount of surface heating – much of which could be re-radiated through an emission line (Thompson, Lyutikov, & Kulkarni 2001).

iv) Torque variations in the SGR/AXP sources. Our picture of torque variations in the SGR and AXP sources is still murky: the level of torque noise is sometimes not much higher than radio pulsars (Kaspi et al. 1999), but at other times the torque can vary by at least a factor $`2`$ over periods of months to years (Woods et al. 2000; Paul et al. 2000; Kaspi et al. 2000). It appears so far that variations in the persistent output (and bursting emission) do not correlate directly with torque variations. Whether there is some hysterisis between the two, or whether they are completely uncorrelated, has important implications for the physical mechanism and can only be determined by redoubling the efforts at phase-connected timing (Kaspi et al. 1999, 2000; Woods et al. 2000). These monitoring campaigns are especially important if, as argued above, the variable torque of the SGRs and AXPs is a distinct physical phenomenon from that encountered in accreting binaries or glitching radio pulsars.

## Acknowledgements

I thank Vicky Kaspi, Shri Kulkarni, and Pete Woods for conversations. Financial support was provided by NASA (NAG 5-3100), the NSERC of Canada, and the Alfred P. Sloan foundation.