text:
# Dark matter in the Galaxy

## 1 Introduction

At some stage of the evolution of the Universe, primordial density fluctuations must have become gravitationally unstable forming dense clumps of dark matter (DM) that have survived until today in the form of galactic halos. In the recent past, galactic halos were successfully modeled as a self-gravitating isothermal gas of particles of arbitrary mass, the density of which scales asymptotically as $`r^2`$, yielding flat rotation curves . The aim of this paper is to describe the halo of our Galaxy in terms of a self-gravitating fermion gas in hydrostatic and thermal equilibrium at finite temperature.

Self-gravitating weakly interacting fermionic matter has been exploited in a wide range of astrophysical phenomena. Originally, self-gravitating degenerate neutrino stars were suggested as a model for quasars , and later neutrino matter was used as a model for dark matter in galactic halos and clusters of galaxies, with a neutrino mass in the $``$ eV range . Recently, degenerate superstars, composed of weakly interacting fermions in the $``$ 10 keV range, were suggested  as an alternative to the supermassive black holes that are believed to exist at the centers of galaxies. It was shown  that such degenerate fermion stars could explain the whole range of supermassive compact dark objects which have been observed so far, with masses ranging from $`10^6`$ to $`3\times 10^9`$$`M_{}`$, merely assuming that a weakly interacting quasistable fermion of mass $`m`$ 15 keV exists in Nature. Most recently, it has been pointed out that a weakly interacting dark matter particle in the mass range $`1\text{ }<\text{ }m/\mathrm{keV}\text{ }<\text{ }5`$ could solve the problem of the excessive structure generated on subgalactic scales in $`N`$-body and hydrodynamical simulations of structure formation in this Universe .

Of course, it is well known that the interval 1-15 keV lies squarely in the cosmologically forbidden mass range for stable active neutrinos $`\nu `$ . The existence of such a massive active neutrino is also disfavored by the Super-Kamiokande data . However, for an initial lepton asymmetry of $`10^3`$, a sterile neutrino $`\nu _s`$ of mass $`m_s10`$ keV may be resonantly produced in the early Universe with near closure density, i.e., $`\mathrm{\Omega }1`$ . The resulting energy spectrum is not thermal but rather cut off so that it approximates a degenerate Fermi gas. In this mass range, sterile neutrinos are also constrained by astrophysical bounds on the radiative $`\nu _s\nu \gamma `$ decay . However, the allowed parameter space includes $`m_s15`$ keV, contributing $`\mathrm{\Omega }_d0.3`$ to the critical density, as favored by the BOOMERANG data . As an alternative possibility, the $``$ 15 keV sterile neutrino could be replaced by the axino  or the gravitino  in soft supersymmetry breaking scenarios.

As the supermassive compact dark objects at the galactic centers are well described by a degenerate gas of fermions, it is tempting to explore the possibility that one could describe both the supermassive compact dark objects and their galactic halos in a unified way in terms of a fermion gas at finite temperature. We will show that this is indeed the case, and that the observed dark matter distribution in the Galactic halo is consistent with the existence of a supermassive compact dark object at the center of the Galaxy which has about the right mass and size, and is in thermal and hydrostatic equilibrium with the halo.

## 2 Standard cosmology

Standard cosmology provides a successful description of the evolution of the Universe from a fraction of a second after the creation until today. A short review on the standard model of cosmology is given in . For our purpose, it is sufficient to state the basic underlying principles. Standard cosmology is based on the following three theoretical assumptions:

1. Cosmological principle. The cosmological principle asserts that the Universe is homogeneous and isotropic on large scales. The most general metric satisfying the cosmological principle is the Friedmann-Robertson-Walker metric

$$ds^2=dt^2a(t)^2\left[\frac{dr^2}{1kr^2}r^2(d\theta ^2+\mathrm{sin}\theta d\varphi ^2)\right],$$
(1)
where the curvature constant $`k`$ takes on the values 1, 0, or -1, for a closed, flat, or open universe, respectively. The time-dependent quantity $`a(t)`$ is the scale factor of the expansion conveniently normalized to unity at present time, i.e., $`a(t_0)=1`$. In other words, $`a`$ is the radius of the Universe measured in units of its current radius.

2. General relativity. Gravity is described by Einstein’s general theory of relativity governed by the equivalence principle and Einstein’s field equations.

3. Perfect fluid. Matter is approximated by a homogeneous perfect fluid. The energy-momentum tensor then takes a simple form

$$T_{\mu \nu }=(\rho +p)u_\mu u_\nu g_{\mu \nu },$$
(2)
where $`\rho `$ and $`p`$ are the density and the pressure of the fluid, respectively.

With these assumptions, the set of Einstein’s equations reduces to the Friedmann-Robertson-Walker (FRW) equations

$$H(t)^2\left(\frac{\dot{a}}{a}\right)^2=\frac{8\pi \rho }{3}\frac{k}{a^2}+\frac{\mathrm{\Lambda }}{3},$$
(3)
$$\frac{\ddot{a}}{a}=\frac{\mathrm{\Lambda }}{3}4\pi (\rho +3p),$$
(4)
where the natural system of units $`\mathrm{}=c=G=1`$ is assumed. The first FRW equation describes the expansion of the Universe. The quantity $`H(t)`$ is the Hubble “constant” and $`\mathrm{\Lambda }`$ is the cosmological constant. We define the critical density as $`\rho _{\mathrm{cr}}=3H_0^2/8\pi `$ and the ratio of the density to the critical density is denoted by $`\mathrm{\Omega }=\rho /\rho _{\mathrm{cr}}`$. The precise present value of the Hubble constant is not known, but the widely accepted value is $`H_0=100h_0\mathrm{kms}^1\mathrm{Mpc}^1`$, with a dimensionless parameter $`h_0`$ between 0.4 and 1. Dividing by $`H(t_0)^2`$, Eq. (3) at $`t=t_0`$ can be conveniently written as a sum rule, i.e.,

$$\mathrm{\Omega }_0k/(H_0a_0)^2+\mathrm{\Omega }_\mathrm{\Lambda }=1,$$
(5)
where $`\mathrm{\Omega }_\mathrm{\Lambda }=\mathrm{\Lambda }/(3H_0^2)`$ is the vacuum energy contribution to the critical density today. Observational evidence favors a flat universe today, i.e., $`k=0`$. Thus, Eq. (5) becomes

$$\mathrm{\Omega }_0+\mathrm{\Omega }_\mathrm{\Lambda }=1.$$
(6)
The second FRW equation (4) describes the acceleration of the expansion. The expansion will accelerate or decelerate, depending on whether the vacuum energy dominates the matter or vice versa. However, even if the cosmological term vanishes, the expansion could accelerate if the dominant component of the DM obeys a peculiar equation of state such that the pressure is negative and higher than one third of the density. One popular example is the scalar field model called quintessence . Another scenario is based on a fluid obeying the Chaplygin gas equation of state $`p1/\rho `$ , which has been intensively investigated for its solubility in 1+1-dimensional space-time, for its supersymmetric extension and connection to $`d`$-branes .

The observational evidence in support of standard cosmology may be summarized in four empirical pillars : Hubble’s law, cosmic background radiation (CBR), anisotropy of CBR, and abundance patterns of light elements. However, despite an overwhelming observational support, a number of problems remain unsolved:

* What caused the Big Bang and the expansion? Did the Universe begin with more than 3+1 dimensions?
* Why is the cosmological constant $`\mathrm{\Lambda }`$ about 50-120 orders of magnitude smaller than the value expected from quantum field theory?
* What caused the initial baryon-antibaryon asymmetry that led to the absence of antimatter today?
* Why is the Universe so smooth on large scales as evidenced by CBR?
* What caused the primordial density fluctuations that provided the seeds for structure formation?
* What does nonbaryonic DM consist of?

A solution to these problems most probably goes beyond the standard model of cosmology and certainly beyond the standard model of particle physics.

## 3 Dark matter

Here, we briefly discuss the DM problem and possible DM candidates. A more detailed analysis may be found in a number of recent review articles .

DM has to be introduced because of the following facts:

* Astronomical observations, such as the flattness of the rotation curves of spiral galaxies and the peculiar motion of galaxies within clusters, strongly indicate that

$$\mathrm{\Omega }_0\frac{\rho _{\mathrm{matt}}}{\rho _{\mathrm{cr}}}\text{ }>\text{ }0.3.$$
(7)
* Consistency with the Big Bang nucleosynthesis implies for baryonic matter

$$0.008h_0^2<\mathrm{\Omega }_B<0.024h_0^2.$$
(8)
* Astronomical observations yield a small relative density of luminous matter

$$\mathrm{\Omega }_{\mathrm{lum}}=0.003h_0^1.$$
(9)

From these facts we conclude that

* About 95% of matter is dark.
* About 60%-90% of DM is nonbaryonic.
* At least 75% of baryonic matter is dark.

Whereas baryonic DM is most likely in the form of relatively standard astrophysical objects, e.g., cold hydrogen clouds or compact objects such as neutron stars, brown dwarfs, MACHOs, and even black holes, the nature of nonbaryonic DM is unknown and still a subject of speculations. Candidates within the standard model are practically excluded and those beyond the standard model have not yet been detected in particle physics experiments. Nevertheless, cosmological and astrophysical observations tell us what these, yet undetected, particles could or could not be.

The different DM scenarios are conveniently classified as hot, warm, and cold DM , depending on the thermal velocities of DM particles in the early Universe.

Hot DM refers to low-mass neutral particles that are still relativistic when galaxy-size masses ($`10^{12}M_{}`$) are first encompassed within the horizon. Hence, fluctuations on galaxy scales are wiped out by the “free streaming” of the dark matter. Standard examples of hot DM are neutrinos and majorons. They are still in thermal equilibrium after the QCD deconfinement transition, which took place at $`T_{\mathrm{QCD}}`$ 150 MeV. Hot DM particles have a cosmological number density comparable with that of microwave background photons, which implies an upper bound to their mass of a few tens of eV.

Warm DM particles are just becoming nonrelativistic when galaxy-size masses enter the horizon. Warm DM particles interact much more weakly than neutrinos. They decouple (i.e., their mean free path first exceeds the horizon size) at $`TT_{\mathrm{QCD}}`$. As a consequence, their number is expected to be roughly an order of magnitude lower and their mass an order of magnitude larger, than hot DM particles. Examples of warm DM are $``$ keV sterile neutrinos, axinos , or gravitinos in soft supersymmetry breaking scenarios . There has been renewed interest in the standard model neutrino as a candidate for warm DM .

Cold DM particles are already nonrelativistic when even globular cluster masses ($`10^6M_{}`$) enter the horizon. Hence, their free streaming is of no cosmological importance. In other words, all cosmologically relevant fluctuations survive in a universe dominated by cold DM. The two main particle candidates for cold dark matter are the lowest supersymmetric weakly interacting massive particles (WIMPs) and the axion.

One of the central issues in dark-matter modeling is the problem of structure formation on subgalactic scales. The combination of cold DM and a small cosmological constant ($`\mathrm{\Lambda }`$CDM) seems to be in good agreement with many observational constraints. However, N-body and hydrodynamical simulations of galaxy formation evidence that $`\mathrm{\Lambda }`$CDM overpredicts structure on small scales . In addition to that, high-resolution simulations generally find a dark matter profile with a central cusp $`\rho r^{1.5}`$ for galactic halos  which seems to contradict the observations.

Clustering on small scales could be suppressed by an upper limit to the phase-space density of DM particles owing either to degeneracy pressure if they are fermions or to a repulsive interaction if they are bosons. The fermion mass would have to lie in the range $`0.1\text{ }<\text{ }m\text{ }<\text{ }10`$ keV. A specific scenario invoking keV mass fermions is the cool-dark-matter proposal  for which candidates exist in shadow-world models and the axino.

## 4 Galactic halo

We now discuss the properties of the halo of our Galaxy assuming that it consists of a self-gravitating gas of keV mass fermions in hydrostatic and thermal equilibrium at finite temperature. The Milky Way Galaxy consists of five major components which are nested within each other . A spheroidal halo with modest concentations of stars and about 170 globular clusters extends out to a radius of perhaps 200 kpc. Within a radius of $`25`$ kpc, the halo contains stars and open clusters that are concentrated into two essentially coplanar disks: the thin disk and the thick disk. In their innermost part the disks merge with a spheroidal bulge, the central concentration of luminous matter in the Galaxy. Finally, a compact dark object with a mass of $`M_\mathrm{c}2.6\times 10^6M_{}`$ is located in the vicinity of the enigmatic radio source Sgr A at the Galactic center , within a radius of 18 mpc.

Here, we demonstrate that by extending the Thomas-Fermi theory to nonzero temperature it is possible to explain, within the same model, both the Galactic halo and the compact dark object at the Galactic center. Extending the Thomas-Fermi theory to finite temperature , it has been shown that, at some critical temperature $`T_\mathrm{c}`$, weakly interacting massive fermionic matter with a total mass below the Oppenheimer-Volkoff limit  undergoes a first-order gravitational phase transition from a diffuse to a clustered state, i.e., a nearly degenerate fermion star. However, during this first-order phase transition a large amount of latent heat must be released in order to substantially decrease the entropy of the initial diffuse configuration. In the absence of a mechanism which would make such a release possible, the system would remain in a thermodynamic quasistable supercooled state close to the point of gravothermal collapse. The Fermi gas will be caught in the supercooled state even if the total mass of the gas exceeds the Oppenheimer-Volkoff limit as a stable condensed state does not exist in this case.

The formation of a quasistable supercooled state may be understood as a process similar to that of violent relaxation, which was introduced to describe rapid virialization of stars of different mass in globular clusters . Through the gravitational collapse of an overdense fluctuation at about one Gyr after the Big Bang, part of gravitational energy transforms into the kinetic energy of random motion of small-scale density fluctuations. The resulting virialized cloud will thus be well approximated by a gravitationally stable thermalized halo. In order to estimate the mass-to-temperature ratio, we assume that an overdense cloud of mass $`M`$ stops expanding at the time $`t_\mathrm{m}`$, reaching its maximal radius $`R_\mathrm{m}`$ and the average density $`\rho _\mathrm{m}=3M/(4\pi R_\mathrm{m}^3)`$. The total energy per particle is just the gravitational energy

$$E=\frac{3}{5}\frac{M}{R_\mathrm{m}}.$$
(10)
From the spherical model of nonlinear collapse  it follows

$$\rho _\mathrm{m}=\frac{9\pi ^2}{16}\overline{\rho }(t_\mathrm{m})=\frac{9\pi ^2}{16}\mathrm{\Omega }_d\rho _{\mathrm{cr}}(1+z_\mathrm{m})^3,$$
(11)
where $`\overline{\rho }(t_\mathrm{m})`$ is the background density at the time $`t_\mathrm{m}`$ or the cosmological redshift $`z_\mathrm{m}`$. We approximate the virialized cloud by a singular isothermal sphere  of the mass of the Galaxy $`M`$ and radius $`R`$. The singular isothermal sphere is characterized by a constant circular velocity $`\mathrm{\Theta }=(2T/m)^{1/2}`$ and the density profile

$$\rho (r)=\frac{\mathrm{\Theta }^2}{4\pi r^2}.$$
(12)
Its total energy per particle is the sum of gravitational and thermal energies, i.e.,

$$E=\frac{1}{4}\frac{M}{R}=\frac{1}{4}\mathrm{\Theta }^2.$$
(13)
Combining Eqs. (10), (11), and (13), we find

$$\mathrm{\Theta }^2=\frac{6\pi }{5}(6\mathrm{\Omega }_d\rho _{\mathrm{cr}}M^2)^{1/3}(1+z_\mathrm{m}).$$
(14)
Taking $`\mathrm{\Omega }_\mathrm{d}=0.3`$, $`M=2\times 10^{12}M_{}`$, $`z_\mathrm{m}=4`$, and $`h_0=0.65`$, we find $`\mathrm{\Theta }220\mathrm{km}\mathrm{s}^1`$, which corresponds to the mass-temperature ratio $`m/T4\times 10^6`$.

Next, we briefly discuss the general-relativistic Thomas-Fermi theory  for a self-gravitating gas, consisting of $`N`$ fermions of mass $`m`$ in equilibrium at a temperature $`T`$, enclosed in a sphere of radius $`R`$. We denote by $`p`$, $`\rho `$, and $`n`$ the pressure, energy density, and particle number density of the gas, respectively. The metric generated by the mass distribution is static, spherically symmetric, and asymptotically flat, i.e.,

$$ds^2=\xi ^2dt^2(12/r)^1dr^2r^2(d\theta ^2+\mathrm{sin}\theta d\varphi ^2).$$
(15)
For numerical convenience, we introduce the parameter

$$\alpha =\frac{\mu }{T}$$
(16)
and the substitution

$$\xi =\frac{\mu }{m}(\phi +1)^{1/2},$$
(17)
where $`\mu `$ is the chemical potential associated with the conserved particle number $`N`$. Using this, the equation of state for a self-gravitating ideal gas may be represented in a parametric form :

$$n=\frac{1}{\pi ^2}_0^{\mathrm{}}𝑑y\frac{y^2}{1+\mathrm{exp}\{[(y^2+1)^{1/2}/(\phi +1)^{1/2}1]\alpha \}},$$
(18)
$$\rho =\frac{1}{\pi ^2}_0^{\mathrm{}}𝑑y\frac{y^2(y^2+1)^{1/2}}{1+\mathrm{exp}\{[(y^2+1)^{1/2}/(\phi +1)^{1/2}1]\alpha \}},$$
(19)
$$p=\frac{1}{3\pi ^2}_0^{\mathrm{}}𝑑y\frac{y^4(y^2+1)^{1/2}}{1+\mathrm{exp}\{[(y^2+1)^{1/2}/(\phi +1)^{1/2}1]\alpha \}}.$$
(20)
We have chosen appropriate length and mass scales $`a`$ and $`b`$, respectively, such that

$$a=b=\sqrt{\frac{2}{g}}\frac{1}{m^2},$$
(21)
or, restoring $`\mathrm{}`$, $`c`$, and $`G`$, we have

$$a=\sqrt{\frac{2}{g}}\frac{\mathrm{}M_{\mathrm{Pl}}}{cm^2}=1.0798\times 10^{10}\sqrt{\frac{2}{g}}\left(\frac{15\mathrm{k}\mathrm{e}\mathrm{V}}{m}\right)^2\mathrm{km},$$
(22)
$$b=\sqrt{\frac{2}{g}}\frac{M_{\mathrm{Pl}}^3}{m^2}=0.7251\times 10^{10}\sqrt{\frac{2}{g}}\left(\frac{15\mathrm{k}\mathrm{e}\mathrm{V}}{m}\right)^2M_{}.$$
(23)
Here, $`M_{\mathrm{Pl}}=\sqrt{\mathrm{}c/G}`$ denotes the Planck mass and $`g`$ the combined spin-degeneracy factor of neutral fermions and antifermions, i.e., g=2 or 4 for Majorana or Dirac fermions, respectively. In this way, the fermion mass, the degeneracy factor, and the chemical potential are eliminated from the equation of state.

Einstein’s field equations for the metric (15) are given by

$$\frac{d\phi }{dr}=2(\phi +1)\frac{+4\pi r^3p}{r(r2)},$$
(24)
$$\frac{d}{dr}=4\pi r^2\rho .$$
(25)
To these two equations we add

$$\frac{d𝒩}{dr}=4\pi r^2(12/r)^{1/2}n,$$
(26)
imposing the particle-number constraint as a condition at the boundary

$$𝒩(R)=N.$$
(27)
Equations (24)-(26) should be integrated using the boundary conditions at the origin:

$$\phi (0)=\phi _0>1;(0)=0;𝒩(0)=0.$$
(28)
It is useful to introduce the degeneracy parameter $`\eta =\alpha \phi /2`$, which, in the Newtonian limit, approaches $`\eta _{\mathrm{nr}}=(\mu _{\mathrm{nr}}V)/T`$. Here, we have introduced the nonrelativistic chemical potential $`\mu _{\mathrm{nr}}=\mu m`$, with $`\mu _{\mathrm{nr}}m`$, and the approximation $`\xi =e^{V/m}1+V/m`$, with $`V`$ being the Newtonian potential. As $`\phi `$ is monotonously decreasing with increasing $`r`$, the strongest degeneracy is obtained at the center with $`\eta _0=\alpha \phi _0/2`$. The parameter $`\eta _0`$, uniquely related to the central density and pressure, will eventually be fixed by the requirement (27). For $`rR`$, the function $`\phi `$ yields the usual empty-space Schwarzschild solution

$$\phi (r)=\frac{\mu ^2}{m^2}\left(1\frac{2M}{r}\right)^11,$$
(29)
with

$$M=(R)=_0^R𝑑r4\pi r^2\rho (r).$$
(30)
Given the temperature $`T`$, the set of self-consistency equations (18)-(26), with the boundary conditions (27)-(30) defines the general-relativistic Thomas-Fermi equation.

The numerical procedure is now straightforward. For a fixed, arbitrarily chosen $`\alpha `$, we first integrate Eqs. (24) and (25) numerically on the interval $`[0,R]`$ and find solutions for various central values $`\eta _0`$. Integrating (26) simultaneously, we obtain $`𝒩(R)`$ as a function of $`\eta _0`$. We then select the value of $`\eta _0`$ for which $`𝒩(R)=N`$. The chemical potential $`\mu `$ corresponding to this particular solution is given by (29). If we now eliminate $`\mu `$ using (16), we finally get the parametric dependence on temperature through $`\alpha `$.

The quantities $`N`$, $`T`$, and, $`R`$ are free parameters in our model and their range and choice are dictated by physics. At $`T=0`$ the number of fermions $`N`$ is restricted by the OV limit $`N_{\mathrm{OV}}=2.89\times 10^9\sqrt{2/g}(15\mathrm{keV}/m)^2M_{}/m`$. However, at nonzero temperature, stable solutions exist with $`N>N_{\mathrm{OV}}`$, depending on temperature and radius. In the following, $`N`$ is required to be of the order $`2\times 10^{12}M_{}/m`$, so that for any $`m`$, the total mass is close to the estimated mass of the halo . As we have demonstrated, the expected particle mass-temperature ratio of the halo is given by $`\alpha m/T=4\times 10^4`$. The halo radius $`R`$ is in principle unlimited; in practice, however, it should not exceed half the average intergalactic distance. It is known that an isothermal configuration has no natural boundary, in contrast to the degenerate case of zero temperature, where for given $`N`$ (up to the OV limit) the radius $`R`$ is naturally fixed by the condition of vanishing pressure and density. At nonzero temperature, with $`R`$ being unbounded, our gas would occupy the entire space, and fixing $`N`$ would make $`p`$ and $`\rho `$ vanish everywhere. Conversely, if we do not fix $`N`$ and integrate the equations on the interval $`[0,\mathrm{})`$, both $`M`$ and $`N`$ will diverge at infinity for $`T>0`$. Thus, one is forced to introduce a cutoff. In an isothermal model of a similar kind , the cutoff was set at the radius $`R`$, where the energy density was by about six orders of magnitude smaller than the central value. Our choice of $`R=200`$ kpc is based on the estimated size of the Galactic halo.

The only remaining free parameters of our model are the fermion mass $`m`$ and the degeneracy factor $`g`$, which always appear in the combination $`m^4g`$. We fix these parameters at $`m=15`$ keV and $`g=2`$, and justify this choice a posteriori.

We now present the results of the calculations for fixed particle number and temperatures near the point of gravothermal collapse. In Fig. 1 the energy per particle defined as $`E=M/Nm`$ is plotted as a function of temperature for fixed $`N=2\times 10^{12}M_{}/m`$. The plot looks very much like that of a canonical Maxwell-Boltzmann ensemble , with one important difference: in the Maxwell-Boltzmann case, the curve continues to spiral inwards ad infinitum approaching the point of the singular isothermal sphere, that is characterized by an infinite central density. In the Fermi-Dirac case, the spiral consists of two, almost identical curves. The inwards winding of the spiral begins for some negative central degeneracy and stops at the point $`T=2.3923\times 10^7m`$, $`E=1.1964\times 10^7m`$, where $`\eta _0`$ becomes zero. This part of the curve, which basically depicts the behavior of a nondegenerate gas, we call the Maxwell-Boltzmann branch. By increasing the central degeneracy parameter further to positive values, the spiral begins to unwind outwards very close to the inwards winding curve. The outwards winding curve will eventually depart from the Maxwell-Boltzmann branch for temperatures $`T\text{ }>\text{ }10^3m`$. Further increase of the central degeneracy parameter brings us to a region where general-relativistic effects become important. The curve will exhibit another spiral for temperatures and energies of the order of a few $`10^3m`$ approaching the limiting temperature $`T_{\mathrm{}}=2.4\times 10^3m`$ and energy $`E_{\mathrm{}}=3.6\times 10^3m`$, with both the central degeneracy parameter and the central density approaching infinite values. It is remarkable that gravitationally stable configurations with arbitrary large central degeneracy parameters exist at finite temperature even though the total mass exceeds the OV limit by several orders of magnitude.

The results of the numerical integration of Eqs. (24) and (25), without restricting $`N`$, are presented in Fig. 2, where we plot the particle number $`N`$ as a function of the central degeneracy parameter $`\eta _0`$ for several values of $`\alpha `$ close to $`4\times 10^6`$. For fixed $`N`$, there is a range of $`\alpha `$, where the Thomas-Fermi equation has multiple solutions. For example, for $`N=2\times 10^{12}`$ and $`\alpha =4\times 10^6`$ six solutions are found, which are denoted by (1), (2), (3), (3’), (2’), and (1’), corresponding to the values $`\eta _0=`$ -30.528, -25.354, -22.390, 29.284, 33.380, and 40.479, respectively.

In Figs. 3 and 4 we plot the corresponding density profiles and enclosed masses, respectively. For negative central value $`\eta _0`$, for which the degeneracy parameter is negative everywhere, the system behaves basically as a Maxwell-Boltzmann isothermal sphere. Positive values of the central degeneracy parameter $`\eta _0`$ are characterized by a pronounced central core of mass of about $`2.5\times 10^6M_{}`$ within a radius of about 20 mpc. The presence of this core is obviously due to the degeneracy pressure of the Fermi-Dirac statistics. A similar structure was obtained in collisionless stellar systems modeled as a nonrelativistic Fermi gas .

Figs. 3 and 4 show two important features. First, a galactic halo at a given temperature $`T`$ may or may not have a central core, depending on whether the central degeneracy parameter $`\eta _0`$ is positive or negative. As the potential is nearly harmonic up to about 1 to 10 kpc for negative $`\eta _0`$, this may favor the formation of a barred galaxy. Second, the closer to zero $`\eta _0`$ is, the smaller the radius is at which the $`r^2`$ asymptotic behavior of the density begins. The flattening of the Galactic rotation curve begins in the range $`1\text{ }<\text{ }r/\mathrm{kpc}\text{ }<\text{ }10`$, hence the solution (3’) most likely describes the Galaxy’s halo. This may be verified by calculating the rotational curves in our model. We know already from our estimate (14) that our model yields the correct asymptotic circular velocity of 220 km/s. In order to make a more realistic comparison with the observed Galactic rotation curve, we must include two additional matter components: the bulge and the disk. The bulge is modeled as a spherically symmetric matter distribution of the form

$$\rho _\mathrm{b}(s)=\frac{e^{hs}}{2s^3}_0^{\mathrm{}}𝑑u\frac{e^{hsu}}{[(u+1)^81]^{1/2}},$$
(31)
where $`s=(r/r_0)^{1/4}`$, $`r_0`$ is the effective radius of the bulge and $`h`$ is a parameter. We adopt $`r_0=2.67`$ kpc and $`h`$ yielding the bulge mass $`M_\mathrm{b}=1.5\times 10^{10}M_{}`$ . In Fig. 5 the mass of halo and bulge enclosed within a given radius is plotted for various $`\eta _0`$. Here, the gravitational backreaction of the bulge on the fermionic halo has been taken into account.

The data points, indicated by squares, are the mass $`M_\mathrm{c}=2.6\times 10^6M_{}`$ within 18 mpc, estimated from the motion of the stars near Sgr A , and the mass $`M_{50}=5.4_{3.6}^{+0.2}\times 10^{11}`$ within 50 kpc, estimated from the motions of satellite galaxies and globular clusters . Variation of the central degeneracy parameter $`\eta _0`$ between 24 and 32 does not change the essential halo features.

In Fig. 6 we plot the circular velocity components of the halo, the bulge, and the disk. The contribution of the disk is modeled as

$$\mathrm{\Theta }_\mathrm{d}(r)^2=\mathrm{\Theta }_\mathrm{d}(r_\mathrm{o})^2\frac{1.97(r/r_\mathrm{o})^{1.22}}{[(r/r_\mathrm{o})^2+0.78^2]^{1.43}},$$
(32)
where we have taken $`r_\mathrm{o}=13.5`$ kpc and $`\mathrm{\Theta }_\mathrm{d}=100`$ km/s. Here it is assumed for simplicity that the disk does not influence the mass distribution of the bulge and the halo. Choosing the central degeneracy $`\eta _0=28`$ for the halo, the data by Merrifield and Olling  are reasonably well fitted.

We now turn to the discussion of our choice of the fermion mass $`m=15`$ keV for the degeneracy factor $`g=2`$. To that end, we investigate how the mass of the central object, i.e., the mass $`M_\mathrm{c}`$ within 18 mpc, depends on $`m`$ in the interval 5 to 25 keV, for various $`\eta _0`$. We find that $`m15`$ keV always gives the maximal value of $`M_\mathrm{c}`$ ranging between 1.7 and 2.3 $`\times 10^6M_{}`$ for $`\eta _0`$ between 20 and 28. Hence, with $`m15`$ keV we get the value closest to the mass of the central object $`M_\mathrm{c}`$ estimated from the motion of the stars near Sgr A .

## 5 Conclusions

In summary, using the Thomas-Fermi theory, we have shown that a weakly interacting fermionic gas at finite temperature yields a mass distribution that successfully describes both the center and the halo of the Galaxy. For a fermion mass $`m15`$ keV, a reasonable fit to the rotation curve is achieved with the temperature $`T=3.75`$ meV and the central degeneracy parameter $`\eta _0=28`$. With the same parameters, we obtain the mass $`M_{50}=5.04\times 10^{11}M_{}`$ and $`M_{200}=2.04\times 10^{12}M_{}`$ within 50 and 200 kpc, respectively. These values agree quite well with the mass estimates based on the motions of satellite galaxies and globular clusters . Moreover, the mass of $`M_\mathrm{c}2.27\times 10^6M_{}`$, enclosed within 18 mpc, agrees reasonably well with the observations of motion of stars near the compact dark object at the center of the Galaxy.

### Acknowledgement

We thank P. Salucci for valuable discussions and comments. We are grateful to M.R. Merrifield and R.P. Olling for sending us the Galactic rotation curve. This research is in part supported by the Foundation of Fundamental Research (FFR) grant number PHY99-01241 and the Research Committee of the University of Cape Town. The work of N.B. is supported in part by the Ministry of Science and Technology of the Republic of Croatia under Contract No. 0098002.


format: text

text:
# Higher Order Hybrid Monte Carlo at Finite Temperature

## 1 Introduction

In lattice QCD the hybrid Monte Carlo (HMC) algorithm  is widely used for simulations of even-flavor of quarks111Odd-flavor simulations of the HMC algorithm are also possible by modifying the Hamiltonian .. These simulations are usually difficult tasks, especially at small quark masses where the computational cost of the matrix solver which is the most time consuming part of the HMC algorithm grows. In order to obtain reliable results within limited computer resources it is important to find an efficient way to implement the HMC algorithm so that the total computational cost of the algorithm is minimized .

The basic idea of the HMC is a combination of 1) molecular dynamics and 2) Metropolis test. Usually the second order leapfrog integrator is used at the molecular dynamics (MD) step. The integrator causes $`𝒪(\mathrm{\Delta }t^3)`$ integration errors, where $`\mathrm{\Delta }t`$ denotes the step-size. Due to the integration errors the Hamiltonian is not conserved. The errors introduced by this integrator have to be removed by the Metropolis test, i.e. accept the new configuration with a probability $`\mathrm{min}(1,\mathrm{exp}(\mathrm{\Delta }H))`$ where $`\mathrm{\Delta }H=H_{new}H_{old}`$ is an energy difference between the starting Hamiltonian $`H_{old}`$ and the new Hamiltonian $`H_{new}`$ at the end of the trajectory.

The acceptance at the Metropolis step depends on the magnitude of the energy difference induced by the integration errors. If a higher order integrator is used at the MD step, the integration errors can be reduced. Therefore one may easily imagine that the performance of the HMC increases with the higher order integrator. However this is not always true since the higher order integrator has more arithmetic operations than the lower one and this might decrease the performance. The total performance should be measured by taking account of two effects: acceptance and the number of arithmetic operations. In Ref. the performance of the higher order integrator at zero temperature ($`\beta =0`$) was studied systematically and it turned out that for the simulation parameters used for the current large-scale simulations with the Wilson fermion action the 2nd order integrator is the best one. The main reason why the higher order integrators are not so effective is that the energy difference caused by the errors of the higher order integrator increases more rapidly than that of the lower one as the quark mass decreases. In Ref. it is shown that at finite temperature the quark mass dependence of the energy difference is small. If so, the conclusion of Ref. might change at finite temperature. In this letter we test higher order integrators at finite temperature and demonstrate that they can actually perform better than the lower order.

## 2 Higher order integrator

In this section we define higher order integrators studied here. Our definition is parallel to that of Ref.. Let $`H`$ be a Hamiltonian given by

$$H=\frac{1}{2}p^2+S(q)$$
(1)
where $`q=(q_1,q_2,\mathrm{})`$ and $`p=(p_1,p_2,\mathrm{})`$ are coordinate variables and conjugate momenta respectively. $`S(q)`$ is a potential term of the system considered. For the lattice QCD, $`q`$ correspond to link variables and $`S(q)`$ consists of gauge and fermion actions.

In the MD step we solve Hamilton’s equations of motion,

$$\{\begin{array}{cc}\frac{dq_i}{dt}\hfill & =\frac{H}{p_i}\hfill \\ \frac{dp_i}{dt}\hfill & =\frac{H}{q_i}.\hfill \end{array}$$
(2)
In general these equations are not solvable analytically. Introducing a step-size $`\mathrm{\Delta }t`$, the discretized version of the equations are solved. In the conventional HMC simulations the 2nd order leapfrog scheme, which causes $`𝒪(\mathrm{\Delta }t^3)`$ step-size error, is used to solve the equations. This scheme is written as

$$\{\begin{array}{cc}q(t+\frac{\mathrm{\Delta }t}{2})\hfill & =q(t)+\frac{\mathrm{\Delta }t}{2}p(t)\hfill \\ p(t+\mathrm{\Delta }t)\hfill & =p(t)\mathrm{\Delta }t\frac{S(q(t+\frac{\mathrm{\Delta }t}{2}))}{q}\hfill \\ q(t+\mathrm{\Delta }t)\hfill & =q(t+\frac{\mathrm{\Delta }t}{2})+\frac{\mathrm{\Delta }t}{2}p(t+\mathrm{\Delta }t).\hfill \end{array}$$
(3)
Eq.(3) forms an elementary MD step. This elementary MD step is performed repeatedly $`N`$ times. The trajectory length $`\tau `$ is given by $`\tau =N\times \mathrm{\Delta }t`$.

Any integrator which satisfies two conditions: 1)time reversible and 2)area preserving can be used for the MD step of the HMC. These conditions are needed to satisfy the detailed balance. When we use the Lie algebraic formalism  we can easily construct higher order integrators which satisfy the above conditions. From the Lie algebraic formalism we find that higher order integrators can be constructed by combining lower order integrators. Let $`G_{2nd}(\mathrm{\Delta }t)`$ be an elementary MD step of the 2nd order integrator with a step-size $`\mathrm{\Delta }t`$. The 4th order integrator $`G_{4th}(\mathrm{\Delta }t)`$ is constructed from a product of $`three`$ 2nd order integrators as

$$G_{4th}(\mathrm{\Delta }t)=G_{2nd}(a_1\mathrm{\Delta }t)G_{2nd}(a_2\mathrm{\Delta }t)G_{2nd}(a_1\mathrm{\Delta }t),$$
(4)
where the coefficients $`a_i`$ are given by

$$a_1=\frac{1}{22^{1/3}},$$
(5)
$$a_2=\frac{2^{1/3}}{22^{1/3}}.$$
(6)
Eq.(4) means that there are three elementary MD steps: (i) first we integrate the equations by eq.(3) with a step-size of $`a_1\mathrm{\Delta }t`$, (ii) then proceed the 2nd order integration with a step-size of $`a_2\mathrm{\Delta }t`$, (iii) finally integrate the equations with a step-size of $`a_1\mathrm{\Delta }t`$. After performing these three elementary MD steps sequentially we obtain the 4th order integrator with the step-size $`\mathrm{\Delta }t`$. This construction scheme can be generalized to the higher $`even`$-order integrators. The $`(2k+2)`$-th order integrator is given recursively by

$$G_{2k+2}(\mathrm{\Delta }t)=G_{2k}(b_1\mathrm{\Delta }t)G_{2k}(b_2\mathrm{\Delta }t)G_{2k}(b_1\mathrm{\Delta }t),$$
(7)
where the coefficient $`b_i`$ are

$$b_1=\frac{1}{22^{1/(2k+1)}}$$
(8)
$$b_2=\frac{2^{1/(2k+1)}}{22^{1/(2k+1)}}.$$
(9)
Compared with the 2nd order integrator, the computational cost of the $`n`$-th order one constructed from eq.(7) grows as $`3^{n/21}`$. For instance the 6th order integrator has 9 times more arithmetic operations than those of the 2nd order one. Yoshida  found $`three`$ parameter sets of 6th order integrators with less arithmetic operations (Table.1). Yoshida’s 6th order integrators consist of $`seven`$ $`G_{2nd}`$’s as

$`G_{6th}(\mathrm{\Delta }t)`$ $`=`$ $`G_{2nd}(w_3\mathrm{\Delta }t)G_{2nd}(w_2\mathrm{\Delta }t)G_{2nd}(w_1\mathrm{\Delta }t)G_{2nd}(w_0\mathrm{\Delta }t)`$ (10)
$`\times G_{2nd}(w_1\mathrm{\Delta }t)G_{2nd}(w_2\mathrm{\Delta }t)G_{2nd}(w_3\mathrm{\Delta }t).`$
In Ref. these 6th order integrators were examined and one of three parameter sets, denoted by $`Y1`$, was found to give smaller integration errors than those of others. In this study the parameter set $`Y1`$ ( see Table.1 ) is used for the 6th order integrator.

## 3 Efficiency of the HMC algorithm

In order to compare among various integrators one needs a criterion which ranks integrators. Following Ref. we utilize the efficiency function $`E_{ff}`$ constructed from a product of acceptance $`P_{acc}`$ and step-size $`\mathrm{\Delta }t`$:

$$E_{ff}=P_{acc}\mathrm{\Delta }t.$$
(11)
This function has one maximum at a certain step-size which we denote $`\mathrm{\Delta }t_{opt}`$. Using the energy difference $`\mathrm{\Delta }H`$ the acceptance is expressed by

$`P_{acc}`$ $`=`$ $`\mathrm{erfc}({\displaystyle \frac{1}{2}}\mathrm{\Delta }H^{1/2}),`$ (12)
where erfc is the complementary error function. In stead of using eq.(12), when $`\mathrm{\Delta }H`$ is small, we may use

$$P_{acc}=\mathrm{exp}(\frac{2}{\sqrt{\pi }}\frac{1}{8}\mathrm{\Delta }H^2^{1/2}).$$
(13)
Although mathematically speaking eq.(13) is valid only for $`\mathrm{\Delta }H^2^{1/2}1`$, the numerical study  shows that eq.(13) approximates the acceptance quit well for $`\mathrm{\Delta }H^2^{1/2}3`$ which corresponds to $`P_{acc}20\%`$. This is enough for our purpose since typically the acceptance of the HMC is taken to be $`P_{acc}50\%`$.

In the lowest order of $`\mathrm{\Delta }t`$, $`\mathrm{\Delta }H^2^{1/2}`$ of the $`n`$-th order integrator is expressed as

$$\mathrm{\Delta }H^2^{1/2}=C_nV^{1/2}\mathrm{\Delta }t^n+𝒪(\mathrm{\Delta }t^{n+1}),$$
(14)
where $`V`$ is volume of the system and $`C_n`$ is a Hamiltonian dependent coefficient.

Substituting eq.(14) into eq.(13) one obtains

$$P_{acc}=\mathrm{exp}(\stackrel{~}{C}_nV^{\frac{1}{2}}\mathrm{\Delta }t^n),$$
(15)
where $`\stackrel{~}{C}_nC_n/\sqrt{2\pi }`$. If one uses eq.(15) for eq.(11) one can easily obtain the optimal step-size which gives the maximum efficiency:

$$\mathrm{\Delta }t_{opt}=\sqrt[n]{\frac{1}{n\stackrel{~}{C}_nV^{\frac{1}{2}}}}.$$
(16)
Furthermore substituting $`\mathrm{\Delta }t_{opt}`$ to eq.(15) one obtains the optimal acceptance222A recent study  shows that if one considers higher order effects of $`\mathrm{\Delta }t`$ the optimal acceptance might be slightly changed. In the current study we omit this small effect. as

$`P_{acc}_{opt}`$ $`=`$ $`\mathrm{exp}({\displaystyle \frac{1}{n}})`$ (17)
$`=`$ $`\{\begin{array}{cc}0.61\hfill & \text{2nd}\hfill \\ 0.78\hfill & \text{4th}\hfill \\ 0.85\hfill & \text{6th}.\hfill \end{array}`$ (21)
Note that the above result does not depend on the specific Hamiltonian and can be applied for any model.

Using eq.(16) and (17) we obtain the optimal efficiency of the $`n`$-th order integrator:

$$(E_{ff})_{opt}^{nth}=\mathrm{exp}(\frac{1}{n})\sqrt[n]{\frac{1}{n\stackrel{~}{C}_nV^{\frac{1}{2}}}},$$
(22)
We use eq.(22) to compare among integrators. Eq.(22) is easy to handle because eq.(22) has one unknown parameter $`\stackrel{~}{C}_n`$ and the value of $`\stackrel{~}{C}_n`$ can be estimated easily from a single simulation on a rather small lattice.

Now let us compare $`n`$-th and $`m`$-th order integrators $`(n>m)`$. If one obtains a gain $`G`$ for the $`n`$-th order integrator over the $`m`$-th order one, the following condition must be satisfied:

$$(E_{ff})_{opt}^{nth}=Gk_{nm}(E_{ff})_{opt}^{mth},$$
(23)
where $`k_{nm}`$ is a relative cost factor needed to implement the $`n`$-th order integrator against the $`m`$-th one and $`G`$ is defined so that both $`n`$-th and $`m`$-th order integrators are equally effective with $`G=1`$. In our case, $`k_{42}=3`$ and $`k_{64}=7/3`$. Substituting eq.(22) to eq.(23) and rewriting the equation one obtains the lattice volume needed to have the gain $`G`$:

$$V^{\frac{n}{2}\frac{m}{2}}=(Gk_{nm}\mathrm{exp}(\frac{1}{m}+\frac{1}{n}))^{nm}\left(\frac{1}{m\stackrel{~}{C}_m}\right)^n(n\stackrel{~}{C}_n)^m.$$
(24)
Our formulas are based on the assumption that $`\mathrm{\Delta }H^2^{1/2}`$ satisfies eq.(14). The validity of eq.(14) depends on the action we take and simulation parameters. In general it is expected that the contribution of the higher order terms to eq.(14) becomes small as the lattice size increases. Let us rewrite eq.(14) as

$$\mathrm{\Delta }H^2^{1/2}=C_nV^{1/2}\mathrm{\Delta }t^n+C_{n+1}V^{1/2}\mathrm{\Delta }t^{n+1},$$
(25)
where $`C_{n+1}V^{1/2}\mathrm{\Delta }t^{n+1}`$ stands for the first relevant higher order term. We measure the contribution of the higher order term by the ratio:

$$\frac{C_{n+1}V^{1/2}\mathrm{\Delta }t^{n+1}}{C_nV^{1/2}\mathrm{\Delta }t^n}=\frac{C_{n+1}}{C_n}\mathrm{\Delta }t.$$
(26)
In order to keep the constant acceptance, $`\mathrm{\Delta }t`$ must be taken to be small as the lattice size increases. Thus eq.(26), ie. the contribution of the higher order term, becomes small for larger lattices.

Fig.1 shows a typical example of $`C_n`$ on a $`16^3\times 4`$ lattice as a function of $`\mathrm{\Delta }t`$. $`\mathrm{\Delta }H^2^{1/2}`$ are well expressed by the functions proportional to $`\mathrm{\Delta }t^n`$ except for large $`\mathrm{\Delta }t`$. We are only interested in the acceptance region indicated by eq.(21): $`P_{acc}60\%`$ which corresponds to $`\mathrm{\Delta }H^2^{1/2}1`$333See Fig.1 in .. It seems that for $`\mathrm{\Delta }H^2^{1/2}1`$ the effect of the higher order terms is small.

## 4 Simulation results

We use the plaquette gauge action and standard Wilson fermion action with two flavors of degenerate quarks ($`n_f=2`$). We first determine the coefficients $`C_n`$ at zero and finite temperature. This can be done by measuring $`\mathrm{\Delta }H^2^{1/2}`$ at a small step-size and substituting the results into eq.(14). The trajectory length $`\tau `$ is set to 1. We choose $`\beta =5.75`$ and make simulations on both $`12^4`$ and $`12^3\times 4`$ lattices. The critical kappa $`\kappa _c`$ of $`n_f=2`$ at this $`\beta `$ is around 0.157444This value is estimated from Figure 5 in Ref.. We use several $`\kappa `$’s in a range of $`0.1\kappa 0.155`$. In this range we maintain the confinement phase on the $`12^4`$ lattice and the deconfinement phase on the $`12^3\times 4`$ lattice. In this study we refer to the results on the $`12^4`$($`12^3\times 4`$) lattice as those at zero(finite) temperature.

Figs.2-4 show $`C_n`$ as a function of $`m_qa`$. Here $`m_qa`$ is defined by $`m_qa=(1/\kappa 1/\kappa _c)/2`$. At large quark masses where the fermionic effects are negligible, values of $`C_n`$ at zero and finite temperature coincide each other. This may indicate that the contribution to $`C_n`$ from the gauge sector is almost same at zero and finite temperature. At small quark masses $`C_n`$ at zero temperature increases more rapidly than those at finite temperature as $`m_qa`$ decreases. The quark mass dependence of $`C_n`$ at finite temperature, compared to that at zero temperature, is small for all the integrators studied here. This behavior is consistent with the result of Ref.. Substituting values of $`C_n`$ at finite temperature into eq.(24) with $`G=1`$, we calculate the lattice size $`L_e`$ ( here $`V=L_e^3\times N_t`$ ). This lattice size $`L_e`$ is the one with which the higher order integrator and the lower order one perform equally. For a lattice size $`L>L_e`$ the higher order integrator is more effective than the lower one. Fig.5 shows $`L_e`$ from comparison between the 2nd and 4th order integrators (2nd vs. 4th) and between the 4th and 6th order ones (4th vs. 6th). For the case of 4th vs. 6th, $`L_e`$ increases as $`m_qa`$ decreases, which we do not appreciate. On the other hand for the 2nd vs. 4th, $`L_e`$ remains less than 20 even at small $`m_qa`$. This result is contrast to that obtained at zero temperature where $`L_e`$ increases as $`m_qa`$ decreases . The above result encourages us to use the 4th order integrator at finite temperature.

The results in Fig.5, however, just show the lattice on which the both integrators are equally effective. To use the higher order integrator in simulations one must obtain some gain over the lower order one. In Fig.6, using eq.(24), we show the expected gain ( the region between the solid lines ) at $`\kappa =0.1525(m_qa0.094)`$ as a function of lattice size $`L`$. To have $`G=2`$ gain ( which means 2 times faster ) a lattice size $`L100`$ is required. This huge lattice size is still not accessible in the current large-scale simulations. Probably the maximum lattice size accessible at the moment is $`L50`$. Therefore we can not expect a large gain from the 4th order integrator even if it is used now. If we use a lattice with $`L50`$, $`G1.5`$ can be achieved. Thus at the level of the current large-scale simulations, we expect to obtain $`G1.5`$.

We also make simulations at $`\kappa =0.1525`$ to confirm that we can actually obtain some gain for the 4th order integrator over the 2nd one. At $`\kappa =0.1525`$, $`L_e`$ is estimated to be 17.5(9). We choose $`18^3\times 4`$ and $`28^3\times 4`$ lattices. On the $`18^3\times 4`$ lattice we expect $`G1`$ and on the $`28^3\times 4`$ lattice, $`G>1`$. The step-size is adjusted so that the acceptance gives a similar value with eq.(21). The gain $`G`$ is calculated by

$$G=\frac{(P_{acc}\times \mathrm{\Delta }t)_{4th}}{3(P_{acc}\times \mathrm{\Delta }t)_{2nd}},$$
(27)
where a factor of 3 in the denominator comes from the relative cost factor $`k_{42}=3`$. Table 2 shows the simulation results. Using these results, we obtain $`G=1.08(3)`$ on the $`18^3\times 4`$ lattice and $`G=1.35(8)`$ on the $`28^3\times 4`$ lattice (see also Fig.6). As expected the gain increases with $`L`$. The result on the $`28^3\times 4`$ lattice is an example showing that the 4th order integrator is more effective than the 2nd order one.

## 5 Conclusions

We have studied higher order integrators for the HMC algorithm with the Wilson fermion action at finite temperature. Contrast to the zero temperature case, the 4th order integrator at finite temperature can be more effective than the 2nd order one. This was demonstrated by the simulations at $`\beta =5.75`$ on $`L^3\times 4`$ lattices. The gain is dependent of the lattice size. It was shown that on the $`28^3\times 4`$ lattice at $`\beta =5.75`$ and $`\kappa =0.1525`$ the 4th order integrator is about $`35\%`$ faster than the 2nd order one.

When large-scale simulations at finite temperature are planned, it is recommended to check which integrator is effective for the lattice considered. This check can be done easily. First we measure $`C_2`$ and $`C_4`$. This first step does not take much computational time since they can be measured on a small lattice. Then substituting values of $`C_2`$ and $`C_4`$ to eq.(24), we obtain a relation between $`G`$ and $`L`$. If we obtain $`G>1`$ on the lattice planned for the simulations, we should consider to use the 4th order integrator.

Eq.(24) is obtained by using the approximation at small $`\mathrm{\Delta }t`$. There might exist the systematic errors caused by the approximation. Although for the present study we considered the errors to be small, for other actions and simulation parameters they might contribute largely. We have to keep in mind that there might exist such contributions depending on the simulation details.

## Acknowledgments

The simulations were done on the NEC SX-5 at INSAM Hiroshima University and at Yukawa Institute. The author would like to thank Atsushi Nakamura for useful discussion and comments. This work was supported by the Grant in Aid for Scientific Research by the Ministry of Education, Culture, Sports, Science and Technology(No.13740164).


text:
# Large Force Fluctuations in a Flowing Granular Medium

## Abstract

We report the characteristics of the temporal fluctuations in the local force delivered to the wall of a 2D hopper by a granular medium flowing through it. The forces are predominantly impulsive at all flow rates for which the flow does not permanently jam. The average impulse delivered to the wall is much larger than the momentum acquired by a single particle under gravity between collisions, reflecting the fact that momentum is transferred to the walls from the bulk of the flow by collisions. At values larger than the average impulse, the probability distribution of impulses is broad and decays exponentially on the scale of the average impulse, just as it does in static granular media. At small impulse values, the probability distribution evolves continuously with flow velocity but does not show a clear signature of the transition from purely collisional flow to intermittently jamming flows. However, the time interval between collisions tends to a power law distribution, $`P(\tau )\tau ^{3/2}`$, thus showing a clear dynamical signature of the approach to jamming.

Sand, rather than a liquid, is chosen for the contents of an hourglass because the rate of efflux from a column of sand flowing under gravity does not depend on the height to which it is filled. Likewise, in a static column of sand, the pressure at a given depth is independent of the height of the column above it Jaeger96 . In both cases, this is because the weight of the sand in the interior of the column is borne by the lateral containing walls. In the case of a static granular column, experiments Dantu&Travers ; Liu1995 ; Howell1999 , theory Coppersmith1996 ; CatesChaos99 , and numerical simulations Liu1995 ; Radjai1996 ; Makse2000  have all shown that stresses in the bulk are transmitted to the walls in a very spatially inhomogeneous fashion. This is demonstrated by two independent observations: first, the distribution of forces, $`P(f)`$, at the boundaries of the medium is exponential, rather than gaussian Mueth1998 ; Miller1996  Secondly, grains that are highly stressed are organized into filamentary, linear structures called ‘force chains’ that carry a large fraction of the stress Dantu&Travers ; Liu1995 ; Howell1999 .

In this article we report experiments directed toward producing a complementary understanding of a flowing granular medium. The force chains in a static medium are unstable to perturbations perpendicular to their length; when jostled by other grains moving in a flow, do they merely become short-lived, or do they melt away entirely? Concomitant to this presumed annealing away of the stress inhomogeneities, does the broad, exponential force distribution in the static case become gaussian when grains move and rearrange? The answers to these questions have broad consequences for the prospects of applying continuum theories to flowing granular matter. If indeed the bounding walls of a medium can communicate with the interior by transient force chains, then any approach to a continuum limit must take into account these long length scales. Previous analyses of the force distributions in a static granular medium Coppersmith1996 , show that an exponential at large forces may be obtained from scalar models in which the weight of a bead is borne by the beads beneath it in some random proportion. While force chains do not naturally emerge from these models it is tempting to speculate that force chains and broad force distributions are related manifestations of the large stress inhomogeneities in a granular medium. This point of view is also suggested by recent simulations that show a simultaneous narrowing of force distributions and a blurring of the force chains Makse2000 , as static granular media are subjected to compression. On the other hand, a recent proposal OHern2001  for a unified description of jamming in thermal as well as non-thermal systems identifies signatures of the approach to a jammed state in the force distribution, $`P(f)`$. In contrast to earlier models, they suggest that loss of mobility is due to the formation of force chains, whose presence is most clearly signalled by the scarcity of unstressed regions (i.e. a dip or plateau in $`P(f)`$ at $`f<\overline{f}`$, the average force) rather than by the exponential tails at large force values. The formation of a plateau in $`P(f)`$ has been also been reported in a recent simulation Silbert2001  of grains flowing down an inclined chute.

In this article we report measurements of the temporal fluctuations of the forces delivered to the side walls of a hopper (shown in Fig.1a) which contains a 2D flow of steel balls with diameter $`d`$=3.125mm. The sides of the hopper are at a fixed angle of $`10^{}`$ to the vertical, so that there are no locations within the cell where the flow stagnates. A transducer (PCB Piezotronics model 209C01) sits flush against one of the sidewalls of the hopper and measures forces normal to that wall. The diameter of the head of the transducer $`d`$; thus there is only one ball against the transducer at any time. For most of the measurements we report, the transducer is located at a constant distance of 5 cm (=16$`d`$)from the opening of the hopper. The flow velocity through the hopper is varied by changing the opening at the bottom of the hopper in several steps from 3$`d`$ to 16$`d`$. Video measurements of $`V_f`$, the average flow velocity at the location of the transducer, were in accord with expected bulk features of the flow. The efflux rate was a constant during the flow at each opening, and showed a power law dependence on the size of the opening, with an exponent of 1.5. Previous experiments in hopper flow have shown a rich variety of dynamics including density waves Baxter1989  and regimes of intermittent flow and shocks Veje1996-97 . These more complex phenomena do not appear in our experimental configuration. Our measurements are similar to those in Ref. 15 however, in those experiments the hopper angle is $`45^{}`$ and the flow is frictional rather than collisional.

The voltage signal from the transducer was sampled at a rate of 100 kHz with 12-bit resolution. Fig 1b shows a section of a typical trace of voltage as a function of time, obtained during a drain. The trace represents the passage of several balls; since the flow is very dense there is no evident signature of individual balls. However, the trace is made up of clearly separated peaks indicating that the forces against the walls are collisional despite the fact that particles visually appear to be in contact At the two lowest flow velocities we report (opening sizes of $`3d`$, and $`3.3d`$), occasional sustained contacts begin to appear in the force trace though the flow does not jam permanently. However, at all other flow velocities, there are three well-separated time scales as shown in Fig. 1b: The average passage time of a ball past the transducer head, $`\tau _p(=d/V_f)`$, is much greater than the average time between collisions, $`\overline{\tau }`$, which in turn is much greater than the width of a peak. Thus a ball makes several distinct collisions with the transducer as it flows past its face.

Peaks in the voltage signal above a noise threshold are identified and calibrated against known impulses to obtain the normal component of the momentum, $`I`$, transferred during the impact. The horizontal dashed line in Fig. 1b shows the average peak height corresponding to an average impulse, $`\overline{I}=8\times 10^6kg.m/s.`$ Since the wall has to provide enough upward momentum to prevent the acceleration of balls along the entire width of the channel, it is not surprising that $`\overline{I}`$ is greater than the vertical momentum that a ball would gain falling in gravity between collisions, which is estimated as $`g\overline{\tau }1\times 10^6kg.m/s`$. In Fig. 1c, we show the variation of the average impulse, $`\overline{I}`$, and the average collision time $`\overline{\tau }`$, as a function of flow velocity, $`V_f`$. As might be expected, $`\overline{I}`$ increases monotonically with flow velocity, $`V_f`$; $`\overline{\tau }`$ also increases modestly with $`V_f`$. In what follows, we discuss the significant features of the distributions of these quantities.

In Fig. 2a we plot the probability distribution of the magnitude of impulses for several flow rates. The distributions, $`P(I)`$, show exponential decay at large $`I`$, just as in the case of static granular media. Furthermore, $`P(I)`$ retains this exponential shape as the opening size is increased from the smallest values that allow sustained flow ($`3d`$) to relatively large openings ($`16d`$). Fig. 2a also includes data taken at two openings with the transducer located further up (10 cm) from the opening where there is much greater crystalline ordering of the beads than at the downstream position where most of the other data are taken. As in static granular media Blair2001 , the exponential force distribution is unchanged by the degree of order in the packing of the beads. Furthermore, the curves in Fig. 2a vary monotonically with the flow velocity at the location of the transducer irrespective of the location of the transducer. When the impulse axis is scaled by the average impulse for each flow velocity, the large forces ($`I>>\overline{I}`$) all collapse and show exponential behaviour as shown in Fig. 2b.

The fact that the broad, exponential distribution of forces seen in the static situation is not narrowed by the flow allows for two possibilities. The first is that some dynamic version of force chains persists in flow, with chains dissolving and forming through the bulk of the flow, just as they do in the quasistatic shearing regime explored in Ref. 4. The second possibility is that force chains and exponential force distributions are unrelated phenomena. The latter possibility is developed in a recent study of the transition to a jammed state where O’Hern et al.OHern2001  show that in simulations of Lennard-Jones potentials as well as in a 2D model of a sheared foam, the distribution of forces at large force is dictated only by the interparticle potential at short distances and is not affected by shear rate or temperature; for steep enough potentials, an approximately exponential form results. Their simulations indicate that the chief signature of the transition to a jammed state (which in the case of granular media is known to have force chains) comes from features in the force distribution at low values of force, in particular, a peak or a plateau develops close to the average value of force as the temperature or shear rate is lowered toward the threshold of jamming. This is a significant new idea, in that it identifies a static quantity - the force distribution - as the sign of a jamming transition. Our data do show some of the features of the scenario of O’Hern et al. As discussed, the large force tail of $`P(I)`$ is exponential at all flow rates. Likewise, as seen in Fig. 2b, the shape of $`P(I)`$ at small values of forces changes with flow rate. However, at slower flow rates, $`P(I)`$ tends to move upwards rather than bend down and form a peak. Since we are only able to reliably resolve impulses above a threshold value dictated by noise and transducer response, we cannot rule out the formation of a peak at extremely low forces, however any such peak can only form at force values well below the average, in contrast to the examples studied by O’Hern et al. This comparison is complicated by the fact that we measure a time-averaged impulse, whereas the simulations of Refs. 11 and 12 measure the instantaneous force. A further caveat in making the comparison to Ref. 12 is that in their geometry, grains interact via enduring, plastic, intergrain contacts just prior to jamming. In our flow geometry, the collisional regime extends almost all the way to the jamming threshold To2001 .

Even though we do not observe a static signature of jamming via $`P(I)`$, we observe a clear dynamical signature of the transition to the jamming distribution in, $`P(\tau )`$, the distribution of the time intervals between collisions. In Figure 3 we show $`P(\tau )`$ plotted against $`\tau `$ on a log-log scale. As shown earlier in Fig. 1c, the measured average time interval $`\overline{\tau }`$ changes very little with $`V_f`$. However, the distribution $`P(\tau )`$, gets steadily broader (Fig. 3b) as the flow velocity is reduced and tends to a power-law $`P(\tau )\tau ^{3/2}`$ (Fig. 3a) at a flow velocity of $`V_f=14.7cm/s`$ corresponding to $`a=3.3d`$. In the experiment, there is potentially a long-time cut-off for this power law set by the transit time of a ball across the transducer face. For the data in Fig. 3b, this time-scale $`20msecs`$, which lies beyond the largest time interval for which we are able to obtain good statistics. However, if the power-law of $`\tau ^{3/2}`$ is indeed the asymptotic shape of the distribution, then the mean time interval $`\overline{\tau }=P(\tau )𝑑\tau `$ tends to diverge just as in a glass transition (even though the average time computed from a finite, albeit large, data set shows a relatively innocuous dependence on $`V_f`$ as in Fig. 1c). At even slower flow velocities (the two curves above the solid line in Fig. 3b), there are increasingly frequent instances of sustained contact: at $`V_f=11.7cm/s`$ and $`14.7cm/s`$ the percentage of the time spent in sustained contact is $`9\%`$ and $`4\%`$, respectively, compared to $`0.3\%`$ at $`V_f=44.5cm/s`$. This leads to ambiguities in defining the time interval between collisions possibly resulting in the dashed curves being even broader than the power-law of $`\tau ^{3/2}`$. A simulation of gravity-driven channel flows Denniston1999  also showed a power-law distribution of collision times with an exponent of close to 3. In that simulation, however, they were able to compute the distribution of collision times for particles in the bulk, averaged over time, whereas in the experiment we only have access to the collisions at a fixed location on the wall and cannot keep track of any correlations that are convected down the flow with individual balls.

In summary, our data show that the collisional regime extends very close to the threshold of jamming. The broad, exponential distribution of forces in the jammed, static case extends right through the jamming transition and up to the largest flow velocities that we are able to explore. The small force end of the force distribution does evolve as the flow is slowed down, but the indication of an imminent jammed state appears via a dynamical quantity, namely, the distribution of collision times, $`P(\tau )`$. The existence of transient force chains in a flowing medium remains unresolved by the measurements reported here, though it is clear that if force chains are intimately connected with exponential force fluctuations, then they must persist in even the most rapid flows we have investigated. We are pursuing this issue further by direct measurement of spatial correlations between different points in the flow.

###### Acknowledgements.

We thank S. Iwanaga, K. Byers, F. Rouyer, S.M. Dragulin for important contributions to the experiment, and S.R. Nagel, A.J. Liu, and O. Narayan for helpful comments. We also gratefully acknowledge support from NSF CAREER DMR-9874833, the UMass MRSEC (DMR-9809365), and the Schultz Foundation.


text:
# X-ray Flashes from Off-axis Gamma-Ray Bursts

## 1 INTRODUCTION

Recently, a new class of X-ray transients has been recognized. The Wide Field Cameras (WFCs) on the BeppoSAX in the X-ray range $`225\mathrm{keV}`$ have detected some Fast X-ray Transients (FXTs) with a duration less than $`10^3\mathrm{s}`$, which are not triggered and not detected by the Gamma Ray Burst Monitor (GRBM) in the $`\gamma `$-ray range $`40700\mathrm{keV}`$ (Heise et al. 2001; see also Strohmayer et al. 1998; Gotthelf, Hamilton & Helfand 1996; Hamilton, Gotthlf & Helfand 1996). In Heise et al. (2001), these FXTs are defined as X-ray flashes (XRFs). This definition of XRFs excludes the X-ray counterparts of the typical Gamma-Ray Bursts (GRBs) including X-ray-rich GRBs. Seventeen XRFs have been observed in the WFCs on the BeppoSAX in about 5 yr, while 49 GRB counterparts have been observed in the same period.

XRFs have the following properties (Heise et al. 2001). (i) The peak flux of the XRFs ranges between $`10^8`$ and $`10^7\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ (Fig. 2 of Heise et al. 2001). The mean peak flux of the XRFs is about a factor of 3 smaller than that of the GRBs. Nine out of 17 XRFs are detected in either the lowest or the lowest two BATSE energy channels ($`2550`$ and $`50100\mathrm{keV}`$; Kippen et al. 2001). (ii) The ratio of the peak flux and the fluence in the X-ray range ($`210\mathrm{keV}`$) and the $`\gamma `$-ray range ($`50300\mathrm{keV}`$) for nine XRFs are shown in Fig. 3 of Heise et al. (2001). The peak flux ratio extends up to a factor of 100, and the fluence ratio extends up to a factor of 20. (iii) The energy spectrum in the range $`225\mathrm{keV}`$ fits with a single power low with the photon index between 1.2 and 3 and the mean of about 2, while the mean photon index of 36 GRBs in the same X-ray band is about 1, with the range between 0.5 and 3. (iv) The duration of the XRFs ranges between $`10\mathrm{s}`$ and $`200\mathrm{s}`$, which is the same order as that of the GRBs. (v) The event rate of the XRFs is estimated as $`100`$ yr-1 since the WFCs observed $`3`$ yr-1 with the covering $`40^{}\times 40^{}`$ (full width to zero response). (vi) The sky distribution is consistent with being isotropic. The spatial distribution is consistent with being homogeneous in Euclidean space since $`V/V_{max}=0.56\pm 0.12`$ (Heise et al. 2000; Schmidt, Higdon & Hueter 1988).

At present, the origin of the XRFs is not known. Heise et al. (2001) have proposed that XRFs could be GRBs at large redshift $`z>5`$, when $`\gamma `$-rays would be shifted into the X-ray range. However, as they have pointed out in their paper, one cannot explain the duration distribution since no time dilation due to cosmological expansion is observed. There is also a possibility that the XRFs could be dirty fireballs or failed GRBs (e.g., Dermer, Chiang & Böttcher1999; Heise et al. 2001; Huang, Dai & Lu 2002).

Ioka & Nakamura (2001) have proposed that the XRFs could be GRBs observed from the large viewing angle as shown in Figure 1 (see also Nakamura 2000). They computed the kinematical dependence of the peak luminosity, the pulse width and the spectral lag of the peak luminosity on the viewing angle $`\theta _v`$ of a jet. For appropriate model parameters they obtained the peak luminosity-spectral lag relation similar to the observed one. They suggested that the viewing angle of the jet might cause various relations in GRBs such as the peak luminosity-variability relation and the luminosity-width relation. Very recently several authors have also suggested that the viewing angle is the key parameter to understand the various properties of the GRBs (Zhang & M$`\stackrel{´}{\mathrm{e}}`$sz$`\stackrel{´}{\mathrm{a}}`$ros 2001; Rossi, Lazzati & Rees 2001; Salmonson & Galama 2001). In this circumstance, it is meaningful to study the off-axis GRB model for the XRFs by Ioka & Nakamura (2001) in more detail.

In this Letter, we will show that the GRBs observed from the large viewing angle possesses the above-mentioned properties (i)-(vi) of the XRFs. In § 2 we describe a simple jet model for the XRFs. In § 3 we consider the peak flux ratio and the fluence ratio (property (ii)). In § 4 we consider the peak flux, the photon index and the event rate (properties (i), (iii) and (v)). § 5 is devoted to discussion (properties (iv) and (vi)).

## 2 EMISSION MODEL OF X-RAY FLASHES

We apply a simple jet model by Ioka & Nakamura (2001) to the XRFs. There are three timescales that determine the temporal pulse structure of the XRFs: the hydrodynamic timescale $`T_{dyn}`$, the cooling timescale $`T_{cool}`$, and the angular spreading timescale $`T_{ang}`$ (Kobayashi, Piran & Sari 1997; Katz 1997; Fenimore, Madras & Nayakshin 1996). Since we consider that XRFs are the GRBs observed from the large viewing angle, we assume $`T_{cool}T_{dyn}T_{ang}`$ as in the case of GRBs (e.g., Piran 1999; Sari, Narayan & Piran 1996). We adopt an instantaneous emission of an infinitesimally thin shell at $`t=t_0`$ and $`r=r_0`$. Then the observed flux of a single pulse at the observed time $`T`$ is given by

$`F_\nu (T)={\displaystyle \frac{2c^2\beta \gamma ^4A_0(r_0/c\beta \gamma ^2)}{D^2}}{\displaystyle \frac{\mathrm{\Delta }\varphi (T)f\left[\nu \gamma (1\beta \mathrm{cos}\theta (T))\right]}{\left[\gamma ^2(1\beta \mathrm{cos}\theta (T))\right]^2}},`$ (1)
where $`1\beta \mathrm{cos}\theta (T)=(c\beta /r_0)(TT_0)`$ and $`T_0=t_0r_0/c\beta `$. The quantity $`A_0`$ determines the normalization of emissivity, and $`f(\nu ^{})`$ represents the spectral shape (for details, see Ioka & Nakamura 2001, Granot, Piran & Sari 1999, and Woods & Loeb 1999). Let the jet opening half-angle and the viewing angle be $`\mathrm{\Delta }\theta `$ and $`\theta _v`$, respectively (see Figure 1). For $`\mathrm{\Delta }\theta >\theta _v\mathrm{and}0<\theta (T)\mathrm{\Delta }\theta \theta _v`$, $`\mathrm{\Delta }\varphi (T)=\pi `$, otherwise $`\mathrm{\Delta }\varphi (T)=\mathrm{cos}^1\left[\frac{\mathrm{cos}\mathrm{\Delta }\theta \mathrm{cos}\theta (T)\mathrm{cos}\theta _v}{\mathrm{sin}\theta _v\mathrm{sin}\theta (T)}\right]`$. For $`\theta _v<\mathrm{\Delta }\theta `$, $`\theta (T)`$ varies from 0 to $`\theta _v+\mathrm{\Delta }\theta `$ while from $`\theta _v\mathrm{\Delta }\theta `$ to $`\theta _v+\mathrm{\Delta }\theta `$ for $`\theta _v>\mathrm{\Delta }\theta `$. In the latter case, $`\mathrm{\Delta }\varphi (T)=0`$ for $`\theta (T)=\theta _v\mathrm{\Delta }\theta `$. A pulse starts at $`T_{start}=T_0+(r_0/c\beta )(1\beta \mathrm{cos}(\mathrm{max}[0,\theta _v\mathrm{\Delta }\theta ]))`$ and ends at $`T_{end}=T_0+(r_0/c\beta )(1\beta \mathrm{cos}(\theta _v+\mathrm{\Delta }\theta ))`$. The spectrum of the GRBs is well approximated by the Band spectrum (Band et al. 1993). In order to have a spectral shape similar to the Band spectrum, we adopt the following form of the spectrum in the comoving frame,

$`f(\nu ^{})=\left({\displaystyle \frac{\nu ^{}}{\nu _0^{}}}\right)^{1+\alpha _B}\left[1+\left({\displaystyle \frac{\nu ^{}}{\nu _0^{}}}\right)^s\right]^{(\beta _B\alpha _B)/s},`$ (2)
where $`\alpha _B`$ ($`\beta _B`$) is the low (high) energy power law index, and $`s`$ describes the smoothness of the transition between the high and low energy. In the GRBs, $`\alpha _B1`$ and $`\beta _B3`$ are typical values (Preece et al. 2000). Equations (1) and (2) are the basic equations to calculate the flux of a single pulse, which depends on 10 parameters for $`\gamma 1`$, $`\theta _v1`$ and $`\mathrm{\Delta }\theta 1`$: $`\gamma \nu _0^{}`$, $`\gamma \theta _v`$, $`\gamma \mathrm{\Delta }\theta `$, $`r_0/c\beta \gamma ^2`$, $`T_0`$, $`\alpha _B`$, $`\beta _B`$, $`s`$, $`D`$, and $`\gamma ^4A_0`$.

In order to study the dependence on the viewing angle $`\theta _v`$, we fix parameters as $`\gamma \mathrm{\Delta }\theta =10`$, $`\alpha _B=1`$, $`\gamma \nu _0^{}=300\mathrm{keV}`$, $`r_0/c\beta \gamma ^2=10\mathrm{s}`$ and $`s=1`$, since typical GRBs have a break energy of $`300\mathrm{keV}`$ (Preece et al. 2000) and a pulse duration of $`10\mathrm{s}`$. Other parameters, i.e., the viewing angle $`\gamma \theta _v`$, the high energy power law index $`\beta _B`$ and the distance $`D`$, are varied depending on circumstances. We fix the amplitude $`\gamma ^4A_0`$ so that the isotropic $`\gamma `$-ray energy $`E_{iso}=4\pi D^2S(202000\mathrm{keV})`$ equals $`10^{53}\mathrm{erg}`$ when $`\beta _B=3.0`$ and $`\gamma \theta _v=0`$. Here $`S(\nu _1\nu _2)=_{T_{start}}^{T_{end}}F(T;\nu _1\nu _2)𝑑T`$ is the fluence in the energy range $`\nu _1\nu _2`$ and $`F(T;\nu _1\nu _2)=_{\nu _1}^{\nu _2}F_\nu (T)𝑑\nu `$ is the flux in the same energy range. The result is

$$A_0=1.2\mathrm{erg}\mathrm{cm}^2\mathrm{Hz}^1\left(\frac{E_{iso}}{10^{53}\mathrm{erg}}\right)\left(\frac{r_0/c\beta \gamma ^2}{10\mathrm{s}}\right)^2\left(\frac{\gamma }{100}\right)^4.$$
(3)
Note that when we adopt $`\gamma =100`$, the opening half-angle of the jet is similar to the observed one, $`\mathrm{\Delta }\theta 0.1`$, and the total energy corrected for geometry is comparable to the observed value, $`(\mathrm{\Delta }\theta )^2E_{iso}10^{51}\mathrm{ergs}`$ (Frail et al. 2001).

## 3 PEAK FLUX RATIO AND FLUENCE RATIO

In this section, we calculate the peak flux ratio $`R_{peak}=F_{peak}(210\mathrm{keV})/F_{peak}(50300\mathrm{keV})`$ and the fluence ratio $`R_{fluence}=S(210\mathrm{keV})/S(50300\mathrm{keV})`$ and compare the results with observations.

Figure 2 shows the peak flux ratio $`R_{peak}`$ and the fluence ratio $`R_{fluence}`$ as a function of the viewing angle $`\gamma \theta _v`$. When the viewing angle $`\theta _v`$ is larger than the opening half-angle $`\mathrm{\Delta }\theta `$, both the peak flux ratio $`R_{peak}`$ and the fluence ratio $`R_{fluence}`$ increase as the viewing angle $`\gamma \theta _v`$ increases. The ratios, $`R_{peak}`$ and $`R_{fluence}`$, increase as the high-energy index $`\beta _B`$ decreases.

We can understand this behavior as follows. As shown in the Appendix, the maximum frequency $`\nu _{max}`$ at which most of the radiation energy is emitted is estimated as $`\nu _{max}\nu _0^{}/\delta `$, where $`\delta \gamma [1\beta \mathrm{cos}(\theta _v\mathrm{\Delta }\theta )][1+\gamma ^2(\theta _v\mathrm{\Delta }\theta )^2]/2\gamma `$ is the Doppler factor and $`\theta _v>\mathrm{\Delta }\theta `$. Thus the maximum frequency $`\nu _{max}`$ decreases as the viewing angle increases. In the following, we consider two observation bands: the lower energy band $`\nu _1\nu _2\mathrm{keV}`$, and the higher energy band $`\nu _3\nu _4\mathrm{keV}`$. The maximum frequency $`\nu _{max}`$ is larger than the highest observed energy $`\nu _4`$ ($`=300\mathrm{keV}`$ in the present case) when $`\gamma \theta _v<\gamma \theta _v^{(4)}\gamma \mathrm{\Delta }\theta +\sqrt{2\gamma \nu _0^{}/\nu _41}`$. In this case, we observe the low energy part of the Band spectrum in equation (2). Since the low energy power law index is $`\alpha _B=1`$, the peak flux ratio $`R_{peak}=F_{peak}(\nu _1\nu _2\mathrm{keV})/F_{peak}(\nu _3\nu _4\mathrm{keV})`$ and the fluence ratio $`R_{fluence}=S(\nu _1\nu _2\mathrm{keV})/S(\nu _3\nu _4\mathrm{keV})`$ are given by $`R_{peak}R_{fluence}(\nu _2/\nu _4)^{2+\alpha _B}`$, where $`\alpha _B>2`$. Similarly, when the maximum frequency $`\nu _{max}`$ is smaller than the lowest observed energy $`\nu _1=2\mathrm{keV}`$, i.e., $`\gamma \theta _v>\gamma \theta _v^{(1)}\gamma \mathrm{\Delta }\theta +\sqrt{2\gamma \nu _0^{}/\nu _11}`$, the peak flux ratio and the fluence ratio are given by $`R_{peak}R_{fluence}(\nu _1/\nu _3)^{2+\beta _B}`$, where $`\beta _B<2`$.

We compare Figure 2 with observations. Observed peak flux ratios extend up to a factor of 100 and observed fluence ratios extend up to a factor of 20 (Fig. 3 of Heise et al. 2001). One can see that when $`\gamma \mathrm{\Delta }\theta =10\gamma \theta _v\gamma \theta _v^{(1)}3\gamma \mathrm{\Delta }\theta `$ and $`4\beta _B2`$, $`R_{peak}`$ and $`R_{fluence}`$ agree with the observational data. Furthermore, Kippen et al. (2002) reported that $`\nu _{max}`$ ranges between about $`2`$ and $`90\mathrm{keV}`$. For our parameters, this can be reproduced if the viewing angle satisfies $`\mathrm{\Delta }\theta \theta _v\theta _v^{(1)}`$.

## 4 PEAK FLUX, PHOTON INDEX AND EVENT RATE

We calculate the peak flux and the photon index in the energy band $`225\mathrm{keV}`$ as a function of the viewing angle $`\gamma \theta _v`$, and plot it in the peak flux $``$ photon index plane. Figure 3 show the results for $`\beta _B=3`$. The distance is varied from $`D=0.01`$ Gpc to $`D=2.1`$ Gpc for our parameters 111 When we consider the effect of cosmology ($`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`h=0.7`$), $`D2\mathrm{Gpc}`$ corresponds to $`z0.4`$. This does not affect our argument qualitatively, but alters the quantitative results up to a factor of 2.. One can see that the photon index increases and the peak flux decreases as the viewing angle $`\gamma \theta _v`$ increases.

As discussed in § 3, we observe the low- (high-) energy part of the Band spectrum in equation (2) when $`\gamma \theta _v<\gamma \theta _v^{(4)}`$ ($`\gamma \theta _v>\gamma \theta _v^{(1)}`$), where $`\nu _4=25\mathrm{keV}`$ and $`\nu _1=2\mathrm{keV}`$. Therefore, the photon index in the energy range $`225\mathrm{keV}`$ is nearly equal to the low- (high-) energy spectral index $`|\alpha _B|=1`$ ($`|\beta _B|=3`$) when $`\gamma \theta _v<\gamma \theta _v^{(4)}14.8`$ ($`\gamma \theta _v>\gamma \theta _v^{(1)}27.3`$). With the analytical estimates in the Appendix, we can also find that the peak flux $`F_{peak}`$ is approximately given by

$`F_{peak}`$ $``$ $`4.3\times 10^6\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2\left({\displaystyle \frac{D}{1\mathrm{Gpc}}}\right)^2\left[1+\gamma ^2(\theta _v\mathrm{\Delta }\theta )^2\right]^{2+\alpha _B}`$ (4)
$`\times `$ $`\left({\displaystyle \frac{r_0/c\beta \gamma ^2}{10\mathrm{s}}}\right)\left({\displaystyle \frac{\gamma \nu _0^{}}{300\mathrm{keV}}}\right)^{1\alpha _B}\left({\displaystyle \frac{\gamma ^4A_0}{1.2\times 10^8\mathrm{erg}\mathrm{s}\mathrm{cm}^2}}\right),`$
when $`\mathrm{\Delta }\theta \theta _v\theta _v^{(4)}`$ (In practice, eq.  can be applied to larger viewing angles $`\gamma \theta _v30`$. We have confirmed that numerical results can be fitted within 5% errors). The peak flux $`F_{peak}`$ is smaller for larger viewing angle. However, if the distances to such sources are small, $`F_{peak}`$ may be comparable to that of the typical GRBs, which have large distances and small viewing angles.

For comparison, we also plot the observed data in the same figures (Fig. 2 of Heise et al. 2001). One can see that the observed XRFs take place within $`2\mathrm{Gpc}`$ and have a viewing angle $`\gamma \mathrm{\Delta }\theta =10\gamma \theta _v\gamma \theta _v^{(1)}3\gamma \mathrm{\Delta }\theta `$.

We roughly estimate the limits in flux sensitivity of the detectors. On the right-hand side of the oblique dashed line, the peak flux in the $`\gamma `$-ray band $`F_{peak}(40700\mathrm{keV})`$ is larger than the limiting sensitivity of the GRBM ($`10^8\mathrm{erg}/\mathrm{s}\mathrm{cm}^2`$), and such events are observed as GRBs, not as XRFs. The vertical dashed line represents the sensitivity of WFCs, where we assume the integration time of $`\mathrm{ms}`$. Therefore, the observed data of XRFs sit in a fairly narrow region surrounded by two dashed lines.

The distance to the farthest XRF $`D_{XRF}`$ gives the observed event rate of the XRFs. The observed event rate $`R_{XRF}`$ can be estimated as $`R_{XRF}=r_{GRB}n_g(4\pi D_{XRF}^3/3)(f_{XRF}/f_{GRB})`$, where $`r_{GRB}`$ and $`n_g`$ are the event rate of the GRBs and the number density of galaxies, respectively. The quantity $`f_{XRF}`$ ($`f_{GRB}`$) is the solid angle subtended by the direction to which the source is observed as the XRF (GRB). From previous discussions, one can find that the emitting thin shell with opening half-angle $`\mathrm{\Delta }\theta `$ is observed as the XRF (GRB) when the viewing angle is within $`\mathrm{\Delta }\theta \theta _v\theta _v^{(1)}3\mathrm{\Delta }\theta `$ ($`0\theta _v\mathrm{\Delta }\theta `$). Therefore the ratio of each solid angle is estimated as $`f_{XRF}/f_{GRB}(3^21^2)/1^2=8`$. Using this value, we obtain

$`R_{XRF}10^2\mathrm{events}\mathrm{yr}^1\left({\displaystyle \frac{r_{GRB}}{5\times 10^8\mathrm{events}\mathrm{yr}^1\mathrm{galaxy}^1}}\right)\left({\displaystyle \frac{D_{XRF}}{2\mathrm{Gpc}}}\right)^3`$
$`\times \left({\displaystyle \frac{n_g}{10^2\mathrm{galaxies}\mathrm{Mpc}^3}}\right)\left({\displaystyle \frac{f_{XRF}/f_{GRB}}{8}}\right),`$ (5)
which is comparable to the observation.

## 5 DISCUSSION

We have shown that the observed data of the XRFs can be reproduced by a simple jet model of the GRBs. This suggests that the XRFs are identical to GRBs. We may say that in the context of our model, nearby GRBs are observed as XRFs when we see them from the off-axis viewing angle. If the distance to the XRFs is much larger than a few Gpc, they cannot be observed since the observed flux is low. This is consistent with the observed value of $`V/V_{max}0.5`$ since the nearby sources distribute homogeneously in Euclidean space.

Our view of the XRFs is different from that of Heise et al. (2001). They have proposed that XRFs could be GRBs at large redshift $`z>5`$, when $`\gamma `$-rays would be shifted into the X-ray range. However, the observed total duration $`T_{90}^{(obs)}`$ cannot be explained. In our model, $`\gamma `$-rays are shifted into the X-ray range by the relativistic beaming effect. The total duration is equal to the lifetime of the central engine and thus does not depend on the viewing angle $`\theta _v`$. Hence the total duration of the XRFs may be similar to that of the GRBs in our model.

We can calculate $`T_{90}`$, the observed duration of a single pulse in the X-ray band ($`225\mathrm{keV}`$). When the viewing angle ranges from $`\gamma \theta _v=10`$ to $`\gamma \theta _v=30`$, the pulse duration is about $`T_{90}303000\mathrm{s}[(r_0/c\beta \gamma ^2)/10\mathrm{s}]`$. This value is comparable but a little bit inconsistent with the observation since the observed pulse duration $`T_{90}`$, which is the order of the angular spreading timescale, should be less than the total duration $`T_{90}^{(obs)}10200\mathrm{sec}`$, which is the time interval between the first and the last emission. This contradiction can be resolved as follows. So far, we have assumed the isotropic energy of the instantaneous emission $`E_{iso}10^{53}\mathrm{erg}`$ and the time unit $`r_0/c\beta \gamma ^210\mathrm{s}`$. The effect of changing the values of these two parameters appears only in the flux normalization $`(\gamma ^4A_0)(r_0/c\beta \gamma ^2)`$ in equation (1). However, one can see that from equation (3), if one rescales these parameters as $`E_{iso}E_{iso}^{}=10^{53}N^1\mathrm{erg}`$ and $`r_0/c\beta \gamma ^2(r_0/c\beta \gamma ^2)^{}=10N^1\mathrm{s}`$, the flux normalization factor is invariant $`(\gamma ^4A_0)(r_0/c\beta \gamma ^2)=[(\gamma ^4A_0)(r_0/c\beta \gamma ^2)]^{}`$, which implies that the result is unchanged. The value of $`N`$ is the number of instantaneous emissions, since we fix the total emission energy as $`E_{iso}^{(tot)}=10^{53}\mathrm{ergs}`$. If we adopt $`N15`$, $`T_{90}`$ of each emission can be less than $`T_{90}^{(obs)}`$.

Ioka & Nakamura (2001) showed that the variability of GRBs is small for large viewing angle. In addition, our model predicts that the number of pulses of the XRFs is smaller than that of typical GRBs. This can be expected from the following discussions. In this Letter, we consider the time-averaged emissions, which means that successive emissions from multiple subjets with the opening half-angle $`\mathrm{\Delta }\theta ^{(j)}\gamma ^1\mathrm{\Delta }\theta /10`$ are approximated by one spontaneous emission caused by a single jet with the viewing angle $`\theta _v`$ and the opening half-angle $`\mathrm{\Delta }\theta `$. Let the viewing angle of each subjet to be $`\theta _v^{(j)}`$. The observed flux (or fluence) in the X-ray band due to the subjets with $`\theta _v^{(j)}\theta _v+\mathrm{\Delta }\theta `$ are much smaller than that with $`\theta _v^{(j)}\theta _v\mathrm{\Delta }\theta `$, and hence negligible. We have confirmed this in the practical calculation. If $`\theta _v\mathrm{\Delta }\theta `$, the emissions of subjets with $`\theta _v^{(j)}\theta _v\mathrm{\Delta }\theta `$ dominates, while if $`\theta _v0`$, in the GRB case, the emissions from almost all subjets may be detected.

We are grateful to the referee for useful comments. This work was supported in part by Grant-in-Aid for Scientific Research of the Japanese Ministry of Education, Culture, Sports, Science and Technology No.00660 (KI), No.11640274 (TN), No.09NP0801 (TN) and No.14047212 (TN).

## Appendix A ANALYTICAL ESTIMATES

(I) $`\nu _{max}`$ and $`(\nu S_\nu )_{max}`$: In equation (1) the typical value of $`\theta (T)`$ is $`(\theta _v\mathrm{\Delta }\theta )`$ when $`\theta _v>\mathrm{\Delta }\theta `$ since the flux peaks soon after the jet edge becomes visible. Since the function $`\nu ^{}f(\nu ^{})`$ in equation (2) takes a maximum at $`\nu _0^{}`$, $`\nu S_\nu `$ takes a maximum at $`\nu _{max}\nu _0^{}/\delta \delta ^1`$ where $`\delta \gamma [1\beta \mathrm{cos}(\theta _v\mathrm{\Delta }\theta )][1+\gamma ^2(\theta _v\mathrm{\Delta }\theta )^2]/2\gamma `$ and $`S_\nu =_{T_{start}}^{T_{end}}F_\nu (T)𝑑T`$. At $`\nu _{max}`$, $`F_\nu `$ in equation (1) is proportional to $`\delta ^2`$ so that we expect $`(\nu S_\nu )_{max}\delta ^3`$ (Ioka & Nakamura 2001). Note here that $`\mathrm{\Delta }\varphi (T)𝑑T`$ depends on $`\theta _v`$ and $`\delta `$ very weakly.

(II) $`T_{ang}`$ and $`\nu F_\nu ^{peak}`$: The pulse duration $`T_{ang}`$ can be estimated as $`T_{ang}(T_{end}T_{start})\theta _v^2\delta `$ for $`\theta _v\mathrm{\Delta }\theta `$, and $`T_{ang}(T_{end}T_{start})\theta _v\delta ^{1/2}`$ for $`\theta _v\mathrm{\Delta }\theta `$. The peak flux $`F_{peak}`$ can be estimated from the relation $`F_{peak}T_{ang}S\delta ^{1+\alpha _B}(\delta ^{1+\beta _B})`$ when the maximum frequency $`\nu _{max}`$ is higher (lower) than the observed frequency.