text:
# Forming Clusters of Galaxies as the Origin of Unidentified GeV Gamma-Ray Sources

## 1. Introduction

The deepest image of the universe in the high energy gamma-ray band beyond 0.1 GeV has been obtained by the EGRET experiment (Hartman et al. 1999). Identified sources of the third EGRET catalog include the Large Magellanic Cloud, five pulsars, and 66 active galactic nuclei (AGNs) of the blazar class. However, over half of the EGRET sources have not yet been identified as known astronomical objects, and their origin is one of the most interesting mysteries in astrophysics. The distribution of these unidentified sources can be interpreted as the sum of the Galactic component along the Galactic disc ($`|b|\stackrel{<}{}40^{}`$) and another isotropic (i.e., likely extragalactic) component (Özel & Thomson 1996; see also Fig. 2 of Mukherjee et al. 1995). Several candidates have been proposed as the origin of the Galactic unidentified sources, including molecular clouds, supernova remnants, massive stars, and radio-quiet pulsars (see e.g., Gehrels & Michelson 1999 and references therein). However, almost no candidate has been proposed to explain the extragalactic unidentified sources, except for undetected AGNs. Recently, Mirabal et al. (2000) have performed comprehensive follow-up observations for one of the high-latitude unidentified EGRET sources (3EG J1835+5918, $`b=25^{}`$) in X-ray, optical, and radio wavebands. They found that any known class of GeV gamma-ray sources including blazars and pulsars cannot be the origin of 3EG J1835+5918, and it suggests that this source belongs to a new class of GeV gamma-ray emitters.

There are 19 unidentified sources with high galactic latitude of $`|b|>45^{}`$ (about 60 in the whole sky), and 7 of them are noted as ‘em’ in the third EGRET catalog, i.e., possibly extended or multiple sources that are inconsistent with single point sources. Although this ‘em’ designation is quite subjective and we should be careful in interpreting this result (see Hartman et al. 1999 for detail), this may indicate that there is an extended and extragalactic population in the unidentified EGRET sources. Recently Gehrels et al. (2000) chose ‘steady’ unidentified sources from the third EGRET catalog, to eliminate any source with high variability typical of flaring AGNs. These steady sources are mostly distributed at low and mid galactic latitude ($`|b|\stackrel{<}{}40^{}`$) and hence should be the Galactic origin, but 7 sources are still located at $`|b|>45^{}`$, that may be steady extragalactic sources whose number is $`24`$ in the whole sky. Therefore, steady astronomical objects with an extended nature are worth being investigated as a possible origin of unidentified EGRET sources.

It is widely believed that the observed structures in the universe have been produced via gravitational instability. Currently the most successful theory of structure formation is the cold dark matter (CDM) scenario, in which the structures grow up hierarchically from small objects into larger ones. When an object collapses gravitationally and virializes, the baryonic matter in the object is heated by shock waves up to the virial temperature. Particles are expected to be accelerated to high energy by shock acceleration, and accelerated electrons scatter the photons of the cosmic microwave background radiation (CMB) to high energy gamma-ray bands by the inverse-Compton mechanism. Existence of such nonthermal electrons is inferred from radio and hard x-ray observations for some clusters of galaxies (e.g., Fusco-Femiano et al. 1999), although the origin of the nonthermal electrons is not yet clear. It has also recently been argued that this radiation process in the intergalactic medium may explain the diffuse extragalactic gamma-ray background radiation (EGRB) observed in the EGRET range (Loeb & Waxman 2000). However, it is still highly speculative and difficult to test whether this process is really the origin of the EGRB, since the contribution by unresolved active galactic nuclei is also of the same order of magnitudes (see, e.g., Mücke & Pohl 2000 and references therein).

On the other hand, if the structure formation is actually an efficient radiation process of gamma-rays, clusters of galaxies should be strong emitters of gamma-rays when they dynamically form, and the detectability of such forming clusters as discrete sources is of great interest as a new probe of structure formation in the universe as well as a test for the scenario proposed by Loeb & Waxman (2000) for EGRB. In this paper we make a theoretical estimate of the number and angular size of such gamma-ray emitting clusters detectable by EGRET, based on the standard theory of structure formation in the CDM universe. We find a few tens of such forming clusters should have already been detected by EGRET. Detectability of such forming clusters in other wavebands such as optical or x-ray bands will be discussed, in comparison with the conventional clusters of galaxies identified in these wavebands. We will also calculate the EGRB spectrum from structure formation and derive a quantitative relation between the higher cut-off photon energy and magnetic field strength.

Throughout this paper, we assume a CDM universe with the density parameter $`\mathrm{\Omega }_0=0.3`$, the cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, the Hubble constant $`h=H_0/(100\text{km/s/Mpc})=0.7`$, the baryon density parameter $`\mathrm{\Omega }_B=0.015h^2`$, and the density fluctuation amplitude $`\sigma _8=1`$. These cosmological parameters are consistent with various observations including those of the CMB fluctuations (e.g. de Bernardis et al. 2000) and the abundance of x-ray clusters of galaxies (e.g. Eke, Cole & Frenk 1996; Kitayama & Suto 1996b).

## 2. Gamma-ray Luminosity and Flux from Forming Clusters

We first estimate the gamma-ray flux of a gravitationally bound object of total mass $`M`$ that virializes at redshift $`z`$. The typical radius $`r_{\mathrm{vir}}`$, density $`\rho _{\mathrm{vir}}`$, circular velocity $`V_c`$, and temperature $`T_{\mathrm{vir}}`$ of the object can be computed from the spherical collapse model (Peebles 1980; Kitayama & Suto 1996b), that is widely used in study of structure formation. The total gravitational energy given to the baryon gas in the forming cluster is given by $`E_{\mathrm{baryon}}(3/4)(\mathrm{\Omega }_B/\mathrm{\Omega }_0)MV_c^2`$. It is reasonable to expect that a fraction $`\xi _e0.05`$ of this energy goes into accelerated electrons, since such a fraction is inferred for acceleration of cosmic ray electrons in a supernova remnant SN 1006 from x-ray and TeV observations (Koyama et al. 1995; Tanimori et al. 1998) and consistent with the energetics among cosmic-rays, turbulent motions, and supernova rate in our Galaxy. It has also been suggested that the diffuse radio and hard x-ray emissions observed in the Coma cluster (and possibly other several clusters) can be attributed to nonthermal electrons with the electron energy fraction of the same order (e.g., Fusco-Femiano et al. 1999). Therefore, we use $`\xi _e=0.05`$ to determine the normalization of electron energy spectrum throughout this paper.

The maximum Lorentz factor of electrons is constrained by the competition of the Fermi acceleration time and cooling time by inverse-Compton (IC) scattering of CMB photons. The acceleration time is given by

$`t_{\mathrm{acc}}`$ $``$ $`{\displaystyle \frac{r_Lc}{V_s^2}}`$ (1)
$`=`$ $`1.6\times 10^4\gamma _eB_{\mu G}^1V_{s,3}^2\mathrm{yr},`$ (2)
where $`r_L=m_e\gamma _e/(eB)`$ is the Larmor radius of electrons, $`\gamma _e`$ the electron Lorentz factor, $`e`$ the electron charge, $`B_{\mu G}=B/(1\mu G)`$ the magnetic field, and $`V_{s,3}=V_s/(10^3\mathrm{km}/\mathrm{s})`$ the shock velocity that is of the same order of magnitudes with the circular velocity of a halo, $`V_c`$. On the other hand, the IC cooling time is

$`t_{IC}`$ $`=`$ $`{\displaystyle \frac{\gamma _em_ec^2}{(4/3)c\sigma _TU_{\mathrm{CMB}}\gamma _e^2}}`$ (3)
$`=`$ $`2.3\times 10^{12}\gamma _e^1(1+z)^4\mathrm{yr},`$ (4)
where $`\sigma _T`$ is the Thomson cross section and $`U_{\mathrm{CMB}}=4.32\times 10^{13}(1+z)^4\mathrm{erg}\mathrm{cm}^3`$ is the CMB energy density. Equating these expressions of $`t_{\mathrm{acc}}`$ and $`t_{\mathrm{IC}}`$, we have the maximum value of $`\gamma _e`$ as $`\gamma _{e,\mathrm{max}}=1.2\times 10^8(1+z)^2B_{\mu G}^{1/2}V_{s,3}`$. We assume the energy distribution of accelerated electrons as a power-law with an exponential cut-off at $`\gamma _{e,\mathrm{max}}`$, i.e., $`dN_e/d\gamma _e\gamma _e^\alpha \mathrm{exp}(\gamma _e/\gamma _{e,\mathrm{max}})`$, with the standard particle acceleration index of $`\alpha 2`$. As mentioned above, the normalization of this spectrum is set by the equation

$$𝑑\gamma _em_ec^2\gamma _e\frac{dN_e}{d\gamma _e}=\xi _eE_{\mathrm{baryon}},$$
(5)
with the parameter $`\xi _e=0.05`$. The observed photon energy $`ϵ_\gamma `$ scattered by electrons is related to $`\gamma _e`$ as $`ϵ_\gamma =(4/3)\gamma _e^2ϵ_{\mathrm{CMB},0}`$, where $`ϵ_{\mathrm{CMB},0}=6.4\times 10^4`$ eV is the mean photon energy of the CMB at $`z=0`$.

The cooling time of electrons corresponding to photon energy $`ϵ_\gamma `$ can be written as $`t_{\mathrm{IC}}=2.1\times 10^6(ϵ_\gamma /\mathrm{GeV})^{1/2}(1+z)^4`$ yr, and this should be compared with the time for the shock wave to propagate the radius of the virialized halo, $`t_{\mathrm{shock}}r_{\mathrm{vir}}/V_s=r_{\mathrm{vir}}/(4V_c/3)`$. Here we have estimated the shock velocity as $`V_s=(4/3)V_c`$, that is a velocity of a strong shock when a material is shocked by a supersonic piston with a velocity of $`V_c`$, i.e., a typical bulk velocity of material in a collapsed halo. By using the spherical collapse model mentioned above to calculate $`r_{\mathrm{vir}}`$ and $`V_c`$ for the halo, this time scale can be written as $`t_{\mathrm{shock}}=(3/4)^{3/2}\pi ^{1/2}(G\rho _{\mathrm{vir}})^{1/2}1.5(1+z)^{3/2}`$ Gyr, that is essentially the dynamical time of the halo. Note that it depends only on the redshift and not on the halo mass. From this argument the cooling time of electrons emitting gamma-rays above 0.1 GeV is always much shorter than the time scale $`t_{\mathrm{shock}}`$ during which the shock is alive and a halo is an active gamma-ray emitter. Hence the total number of gamma-rays emitted from a forming halo during the time $`t_{\mathrm{shock}}`$ is given as:

$`{\displaystyle \frac{dN_\gamma (ϵ_\gamma ;M,z)}{dϵ_\gamma }}={\displaystyle \frac{m_e\gamma _e}{ϵ_\gamma }}{\displaystyle \frac{dN_e}{d\gamma _e}}{\displaystyle \frac{d\gamma _e}{dϵ_\gamma }},`$ (6)
and the observed photon flux of gamma-rays during the shock propagation time is

$`{\displaystyle \frac{dF(ϵ_\gamma ;M,z)}{dϵ_\gamma }}={\displaystyle \frac{(1+z)}{4\pi d_L^2}}{\displaystyle \frac{dN_\gamma }{dϵ_\gamma }}{\displaystyle \frac{1}{t_{\mathrm{shock}}}},`$ (7)
where $`d_L`$ is the standard luminosity distance.

We introduce a parameter $`\xi _B`$ to determine the magnetic field, which is the fraction of magnetic energy density in the total baryon energy density of the halo, i.e., $`B^2/(8\pi )=\xi _BE_{\mathrm{baryon}}/(4\pi r_{\mathrm{vir}}^3/3)`$. Then by using the spherical collapse model again, the magnetic field of a cluster can be written as

$`B0.17\left({\displaystyle \frac{\xi _B}{10^3}}\right)^{1/2}\left({\displaystyle \frac{M}{10^{15}M_{}}}\right)^{1/3}(1+z)^2\mu \mathrm{G}.`$ (8)
A magnetic field of $``$ 0.1–1 $`\mu `$G is often observed in intracluster medium of rich clusters (Kronberg 1994; Fusco-Femiano et al. 1999; Rephaeli, Gruber, & Blanco 1999), and hence we use $`\xi _B=10^3`$ to be consistent with the observations. It should be noted, however, that this parameter is important only for the maximum photon energy of the gamma-ray spectrum (well beyond 10–100 GeV with $`\xi _B=10^3`$, see §4), and the gamma-ray flux above 100 MeV measured by the EGRET is almost insensitive to this uncertain parameter, unless the particle index $`\alpha `$ significantly deviates from the standard value of 2.

It is interesting to apply the above model to nearby known clusters of galaxies. For example, the Coma cluster has the total mass of $`M10^{15}M_{}`$ and is located at $`z=0.023`$. The flux with these parameters becomes $`F(>0.1\mathrm{GeV})6.5\times 10^7\mathrm{photons}\mathrm{cm}^2\mathrm{sec}^1`$, that is about 10 times brighter than the observational upper limit on this cluster by EGRET, $`4\times 10^8\mathrm{photons}\mathrm{cm}^2\mathrm{sec}^1`$ (Sreekumar et al. 1996). This does not mean, however, that our model is incorrect. Our model is relevant for just-forming clusters of galaxies in which the violent shock generated by gravitational collapse is still alive. The Coma cluster is thought to have formed more than a few dynamical times ago and is now hydrodynamically stable after violent shock has disappeared. Gamma-ray flux can then be much weaker than our estimation. On the other hand, it suggests that there was an epoch during which this cluster was a strong gamma-ray emitter, and also that there may be other clusters visible by EGRET, that are just dynamically forming and have not yet reached hydrostatic equilibrium.

## 3. Expected Number of Gamma-Ray Clusters Detectable by EGRET

The number of such forming clusters of galaxies with flux stronger than $`F`$ can be calculated as

$`N(>F)={\displaystyle 𝑑z_{M(z;F)}^{\mathrm{}}𝑑M\frac{dV}{dz}R_{\mathrm{form}}(M,z)t_{\mathrm{shock}}},`$ (9)
where $`dV/dz`$ is the comoving volume element of the universe, $`R_{\mathrm{form}}`$ the formation rate of dark haloes (or clusters) per unit mass, cosmic time, and comoving volume, and $`M(z;F)`$ the mass of a cluster collapsing at redshift $`z`$ whose flux is $`F`$. Here we have taken into account that clusters are active gamma-ray emitters only during the time $`t_{\mathrm{shock}}`$.

There are several formulae to calculate the formation rate $`R_{\mathrm{form}}`$ in the framework of the standard theory of structure formation. As is well known, the Press & Schechter (1974, PS; Peebles 1980) formalism provides a formula of mass function (i.e., number density of haloes as function of mass and redshift), that is in reasonable agreement with $`N`$-body simulations. Here we want the formation rate of haloes rather than the mass function at a given epoch, because we need to calculate the number of collapsing objects experiencing shock at each epoch. A naive prescription to obtain this quantity is to take a time derivative of the PS mass function, $`R_{dPS}`$, although this is not exactly $`R_{\mathrm{form}}`$ but rather interpreted as $`R_{dPS}=R_{\mathrm{form}}R_{\mathrm{dest}}`$, where $`R_{\mathrm{dest}}`$ is the rate of destruction of haloes by merging into even larger structures. Consequently, $`R_{\mathrm{dPS}}`$ becomes negative at small mass scales where $`R_{\mathrm{dest}}`$ is significant. As shown below, however, the number of objects visible by EGRET is dominated by massive clusters forming in recent past $`(z1)`$ that are the largest structures in the universe, and hence $`R_{\mathrm{dest}}`$ is negligible. Therefore it is a reasonable approximation to use $`R_{dPS}`$ in the mass range where it is positive. Alternatively, there are several extensions to the PS theory for computing $`R_{\mathrm{form}}`$ (e.g., Lacey & Cole 1993; Sasaki 1994; Kitayama & Suto 1996a; Percival & Miller 1999). In what follows we use $`R_{dPS}`$ and the formulae of Sasaki (1994) and Kitayama & Suto (1996a) to take account of theoretical uncertainties in $`R_{\mathrm{form}}`$.

Figure 1 shows the theoretically predicted $`\mathrm{log}N`$-$`\mathrm{log}F`$ of forming clusters. At least a few tens of clusters should be visible by the EGRET, and a significant fraction of the isotropic unidentified EGRET sources can be accounted for. It is also interesting to note that the predicted number is similar to that of ‘em’ isotropic unidentified EGRET sources, i.e., possibly extended sources. The predicted number at the EGRET flux limit also agrees with the number of ‘steady’ unidentified sources with $`|b|>45^{}`$ defined by Gehrels et al. (2000). This result is robust against changes in the adopted prescription for $`R_{\mathrm{form}}`$. In the lower panel of Fig. 1, we also show mean mass, redshift, and angular radius $`\theta _{\mathrm{vir}}`$ corresponding to $`r_{\mathrm{vir}}`$ of such gamma-ray clusters brighter than a given flux. These quantities are $`M10^{15}M_{}`$, $`z0.05`$, and $`\theta _{\mathrm{vir}}1^{}`$ for clusters above the EGRET sensitivity limit. Considering the EGRET angular resolution, the typical radius of $`1^{}`$ is consistent with the fact that a significant fraction of isotropic unidentified sources are indicated as possibly extended.

It is predicted that more than a few thousands of forming clusters will be detected by future missions such as the GLAST (Gehrels & Michelson 1999), and the improved angular resolution may reveal the extended profile for nearby gamma-ray clusters with higher statistical significance. Another important prediction is that GLAST will observe the flattening of the $`\mathrm{log}N`$-$`\mathrm{log}F`$ curve due to the cosmological effects, compared with the expectation of a uniform source distribution in the Euclidean space (dotted line in the upper panel of Fig. 1).

## 4. Extragalactic Gamma-ray Background

Our formulation also allows us to calculate the EGRB flux and spectrum as

$`{\displaystyle \frac{dn_\gamma }{dϵ_\gamma }}={\displaystyle 𝑑z𝑑M\frac{dN_\gamma (ϵ_\gamma ;M,z)}{dϵ_\gamma }R_{\mathrm{form}}(M,z)\frac{dt}{dz}},`$ (10)
where $`(dn_\gamma /dϵ_\gamma )`$ is the gamma-ray number density that is related to the EGRB flux as $`(dF/dϵ_\gamma )=c(dn_\gamma /dϵ_\gamma )/(4\pi )`$, and $`t`$ is the cosmic time. This flux becomes $`ϵ_\gamma ^2(dF/dϵ_\gamma )1\mathrm{k}\mathrm{e}\mathrm{V}\mathrm{cm}^2\mathrm{sec}^1\mathrm{sr}^1`$ at 100 MeV, in good agreement with the observation (Sreekumar et al. 1998) as well as the previous simpler estimation assuming that the average temperature of baryons in the universe is $``$ keV at present (Loeb & Waxman 2000). In fact, we have checked that the mass-averaged temperature of virialized haloes in the universe as a function redshift, that is calculated by the PS theory, agrees within a factor of 2 with a numerical simulation (Cen & Ostriker 1999) on which the previous EGRB estimate was based. This fact gives a justification for the use of the PS theory to calculate the gamma-ray emitting objects. The dotted line in the lower panel of Fig. 1 shows the contribution to the EGRB at 100 MeV by objects brighter than a given flux. We predict that GLAST will resolve about 20–30 % of the EGRB as discrete gamma-ray clusters, if structure formation is the major origin of the EGRB.

The strength of magnetic field in the shocked baryons is important for the question whether the EGRB spectrum extends up to $``$ 100 GeV as observed. In Fig. 2 we show the EGRB spectrum with several values of $`\xi _B`$. This result shows that the magnetic field strength corresponding to $`\xi _B10^5`$ of the baryon energy density is sufficient for the EGRB to extend beyond 100 GeV. The magnetic field observed in the intracluster gas ($`\xi _B10^3`$) is much stronger than this, and it is also theoretically reasonable to expect that the turbulent motion in collapsed objects amplifies the seed magnetic field made by the battery mechanism well beyond $`\xi _B\stackrel{>}{}10^5`$ within the dynamical time (Kulsrud et al. 1997). Therefore, physically reasonable magnetic field strength can explain the extension of the EGRB spectrum beyond 100 GeV, and it is likely that a considerable fraction of gamma-rays above 100 GeV is absorbed by the interaction with the cosmic infrared background radiation producing electron-positron pairs. The effect of intergalactic absorption is shown in the lower panel of Fig 2, using the optical depth of intergalactic pair-production presented in Totani (2000). The absorbed TeV gamma-rays will be reprocessed into GeV gamma-rays by the produced pairs, and distort the EGRB spectrum (e.g., Coppi & Aharonian 1997). Although these secondary GeV gamma-rays are not taken into account here and it is beyond the scope of this paper, it is important to study how large is this spectral distortion in future work.

## 5. Discussion

We here discuss the expected properties of gamma-ray clusters of galaxies. Perhaps the most natural question in this regard would be “Are they already observed in other wavebands such as x-rays or optical surveys?” We have checked that there is no statistically significant association of the ROSAT Brightest Cluster Sample (RBCS, Ebeling et al. 1998), within the 95% error circles of the unidentified sources with $`|b|>30^{}`$ in the EGRET catalog. We have also checked the correlation with the clusters in the revised Abell catalog (Abell, Corwin, & Olowin 1989), and no statistically significant associations are found, either. However, in the following we argue that the gamma-ray clusters proposed in this paper are very difficult to detect in x-rays or optical bands compared with ordinary clusters identified in these wavebands, and hence our scenario is not rejected by these results.

### 5.1. Detectability of gamma-ray clusters in X-rays

We first estimate the expected x-ray flux from gamma-ray clusters. Baryonic gas in most clusters of galaxies observed in x-rays seems to be in approximate hydrostatic equilibrium with the surface brightness well fitted by a density profile, $`\rho _{\mathrm{gas}}(r)[1+(r/r_c)^2]^1`$ (e.g., Sarazin 1988), where $`r_c`$ is the core radius that is typically about $`10`$ times smaller than the virial radius. Since the x-ray emissivity is proportional to $`\rho _{\mathrm{gas}}^2`$, the x-ray emission is strongly concentrated into the central region. Assuming the above density profile and the self-similar model as described in Kitayama & Suto (1996b), <sup>1</sup><sup>1</sup>1Here we have assumed that the core radius is proportional to the virial radius as $`r_c=0.21(h/0.7)^1\mathrm{Mpc}\times [r_{\mathrm{vir}}(M,z)/r_{\mathrm{vir}}(10^{15}M_{},0)]`$ with $`r_{\mathrm{vir}}(10^{15}M_{},0)`$ = 2.6Mpc, where the normalization is chosen to match the local observations (Abramopoulos & Ku 1983; Jones & Forman 1984). a typical gamma-ray cluster detectable by the EGRET with $`M10^{15}M_{}`$ and $`z0.05`$ would have the x-ray flux $`2.4\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in 0.1–2.4 keV.

The inverse-Compton flux is also expected to be comparable with the thermal emission. By equating $`t_{IC}`$ and $`t_{\mathrm{shock}}`$ in §2, we get the cooling photon energy $`ϵ_{\gamma ,\mathrm{cool}}=2.0(1+z)^5`$ keV, below which the electron cooling time is longer than the dynamical time. Then the IC spectrum extends down to around x-ray band with $`dN_\gamma /dϵ_\gamma ϵ_\gamma ^2`$, while it becomes harder at wavelengths longer than x-rays with $`dN_\gamma /dϵ_\gamma ϵ_\gamma ^{1.5}`$. If the gamma-ray flux at 100 MeV is $`10^7\mathrm{photons}\mathrm{cm}^2\mathrm{s}^1`$ that is the EGRET threshold, the IC x-ray flux ($`\nu F_\nu `$) is $`1.6\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. Therefore, the thermal and IC fluxes are well above the flux limit $`4\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ of the RBCS.

However, it takes nearly the dynamical time for the cluster gas to reach hydrostatic equilibrium after the collapse, and gamma-rays from the shock generated by the gravitational collapse are radiated away within that period. Then it is likely that the density profile of gamma-ray emitting clusters is more irregular and extended than ordinary x-ray clusters. In fact, if the unidentified ‘em’ sources in the EGRET catalog are actually extended, they must have typical angular size of about degree, from the source location accuracy of the EGRET. As we have shown, angular size of about 1 is theoretically reasonable if the emission is extended to the virial radius. When the density profile is not concentrated into the central region but rather constant within the virial radius, the x-ray luminosity becomes lower than the self-similar model by a factor of $`3.7`$ because of the lower central density.

Furthermore, the surface brightness of such loose clusters should be drastically dimmer than ordinary x-ray clusters. In the self-similar model with $`r_{\mathrm{vir}}r_c`$, the core gas density is $`\rho _{\mathrm{gas},\mathrm{c}}(1/3)(r_{\mathrm{vir}}/r_c)^2\rho _{\mathrm{gas},\mathrm{vir}}50\rho _{\mathrm{gas},\mathrm{vir}}`$, where $`\rho _{\mathrm{gas},\mathrm{vir}}=(\mathrm{\Omega }_B/\mathrm{\Omega }_0)\rho _{\mathrm{vir}}`$ is the virial gas density that is the average gas density within $`r_{\mathrm{vir}}`$. On the other hand, if the gas density profile of gamma-ray clusters is roughly constant at $`\rho _{\mathrm{gas},\mathrm{vir}}`$ out to $`r_{\mathrm{vir}}`$, the x-ray surface brightness of such a loose cluster is dimmer than the central surface brightness of the self-similar model by a factor of $`(r_c/r_{\mathrm{vir}})(\rho _{\mathrm{gas},\mathrm{c}}/\rho _{\mathrm{gas},\mathrm{vir}})^2200`$, since the x-ray emissivity is proportional to $`\rho _{\mathrm{gas}}^2`$. It crucially affects the detectability of x-rays from gamma-ray clusters. The detectability of x-rays should be described by the signal-to-noise ratio $`(S/N)`$ against the x-ray background flux. The noise level is proportional to (image area)<sup>1/2</sup>, and hence $`S/NF/r`$, where $`F`$ and $`r`$ are the flux and the image radius, respectively. We have compared the value of $`F/r`$ of the extended gamma-ray clusters detectable by the EGRET and those of the clusters in the RBCS. We found that the $`F/r`$ of gamma-ray clusters is by a factor of 3 smaller than the minimum $`F/r`$ of the RBCS clusters. The absence of association between the RBCS and the EGRET sources is therefore not in contradiction to our scenario. On the other hand, deeper observation of candidate gamma-ray clusters by Newton, for example, might detect the x-ray emission extended to about $`1^{}`$ with the flux estimated above, that would provide a clear test of our scenario. Such x-ray emission should reflect the structure of shocks in dynamically forming clusters, and imaging study is of great interest.

### 5.2. Detectability in the Optical Surveys

Here we again emphasize that the gamma-ray clusters are expected to be more extended than clusters that have already stabilized. It is known that the surface density profile of galaxies in a cluster can well be described by the King profile, $`\sigma (r)[1+(r/r_c)^2]^1`$ with the core radius of $``$ 100 kpc that is comparable with the core radius of x-ray profile (e.g., Adami et al. 1998). If we assume a roughly constant surface density out to $`r_{\mathrm{vir}}`$ rather than the King profile for gamma-ray clusters, the average surface density $`\sigma _{\mathrm{av}}N_{\mathrm{gal}}/(\pi r_{\mathrm{vir}}^2)`$ is lower than the central surface density of the King profile $`\sigma _cN_{\mathrm{gal}}/[2\pi r_c^2\mathrm{ln}(r_{\mathrm{vir}}/r_c)]`$, by a factor of $`\sigma _c/\sigma _{\mathrm{av}}[2\mathrm{ln}(r_{\mathrm{vir}}/r_c)]^1(r_{\mathrm{vir}}/r_c)^230`$. Here $`N_{\mathrm{gal}}`$ is the total number of galaxies within $`r_{\mathrm{vir}}`$. This dimming factor is not so significant as that for the x-ray surface brightness, but that should make the optical identification very difficult because of the contamination by foreground and/or background field galaxies. Therefore, we consider that no statistically significant association with the known optically identified clusters does not immediately reject our scenario. Instead, it is necessary to study in the future the correlation between the EGRET sources and galaxy catalogs taking into account the possibility that the gamma-ray clusters are considerably extended. Search in optical bands has an advantage over the search in x-rays, in a sense that the dimming of surface number density compared with ordinary clusters is less severe than x-rays whose emissivity is proportional to $`\rho _{\mathrm{gas}}^2`$. The typical density of such loose clusters is close to the virial density, that is about a few hundreds times higher than the mean density of the universe.

We have also noticed that there are a considerable number of ‘em’ sources in the EGRET sources identified as AGNs. If they were actually extended sources, it might be speculated that some of them are also gamma-ray clusters including an AGN as a member galaxy. Time variability of these sources would be an important test to check this possibility.

### 5.3. On the recent follow-up observations for 3EG J1835+5918

Recent follow-up observations by Mirabal et al. (2000) for one of the high-latitude unidentified EGRET sources (3EG J1835+5918) have found a diffuse x-ray emission from an uncatalogued cluster of galaxies at $`z=0.102`$. Although this cluster is outside the 99% error ellipse of 3EG J1835+5918 whose radius is 12, the separation between the centers of the x-ray cluster and 3EG J1835+5918 is about 0.65, that is within our expectation of the typical angular radius of gamma-ray clusters detectable by the EGRET, $`1^{}`$. As discussed above, x-ray emission is expected from a region where the intracluster gas reached hydrodynamical equillibrium, while gamma-rays are emitted from a region still hydrodynamically unstable. Therefore, it is not surprising that the positions of x-ray and gamma-ray emissions are different unless the separation is well beyond the virial radius of $`1^{}`$.

3EG J1835+5918 is not an ‘em’ source, and this source may not be an extended source. The radius of the 99% confidence ellipse, $`0.2^{}`$, is considerably smaller than the expected angular size of gamma-ray clusters. However, high energy electrons emitting GeV gamma-rays have very short life time compared with the shock propagation time ($`t_{\mathrm{IC}}10^3t_{\mathrm{shock}}`$), and gamma-ray emitting region may be very clumpy in a cluster. (On the other hand, x-ray emitting electrons have a cooling time comparable with or longer than $`t_{\mathrm{shock}}`$ for IC and thermal radiations, and hence x-ray emitting region should be much less clumpy and extended with the size $`\theta _{\mathrm{vir}}`$, as discussed in §5.1.) Therefore, it is possible that the gamma-ray size of the 3EG J1835+5918 is considerably smaller than the physical size of a whole forming cluster. This consideration also suggests a possibility that some of gamma-ray clusters may be observed as multiple sources within $`\theta _{\mathrm{vir}}`$, that may be revealed by future gamma-ray missions.

One of the characteristics that the source of 3EG J1835+5918 must have is very weak radio flux that is at least two orders of magnitudes fainter than any of the securely identified EGRET blazars (Mirabal et al. 2000). The spectrum of blazars is well understood by the two components of radiation by the same population of nonthermal electrons, i.e., synchrotron radiation in radio, optical, and X-ray bands, and inverse-Compton radiation in GeV and TeV gamma-ray bands (e.g., Inoue & Takahara 1996; Kataoka et al. 1999). The ratio of luminosities by the two processes is, as is well-known, given by the ratio of magnetic energy density to the target photon energy density, that is typically of order unity for blazars. On the other hand, this ratio is $`U_B/U_{\mathrm{CMB}}2.7\times 10^3(\xi _B/10^3)(M/10^{15}M_{})^{2/3}`$ for gamma-ray clusters, by using the expression of $`B`$ given in §2. Therefore, the $`U_B/U_{\mathrm{CMB}}`$ ratio is generally much smaller than the unity, and very weak radio flux compared with identified blazars can be reasonably explained.

Based on the above arguments, we suggest that the uncatalogued x-ray cluster near 3EG J1835+5918 may be a gamma-ray cluster proposed in this paper, and further observations for this cluster and the surrounding region are very important. Our model predicts that a cluster emitting a flux of $`6.06\times 10^7\mathrm{photons}\mathrm{cm}^2\mathrm{s}^1`$ above 100 MeV (the flux of 3EG J1835+5918, Hartman et al. 1999) at $`z`$ = 0.102 should have a total mass of $`7\times 10^{15}M_{}`$ and $`r_{\mathrm{vir}}5`$ Mpc ($`\theta _{\mathrm{vir}}0.7^{}`$).

## 6. Conclusions

In this paper we have proposed a new candidate of unidentified EGRET sources: gamma-ray clusters that are just dynamically forming and emit gamma-rays due to inverse-Compton scattering of CMB photons by shock-accelerated electrons. Based on the standard theory of structure formation and assuming the injection of $``$5% of shock energy at the formation into nonthermal electrons, we have shown that a few tens of such clusters should have already been detected by EGERT, and a significant fraction of the isotropic component of unidentified EGRET sources can be accounted for. Such gamma-ray clusters are expected to be very extended; the x-ray surface brightness and surface number density of galaxies could be lower than those of ordinary clusters by a factor of $`200`$ and $`30`$, respectively. Therefore it should have been very difficult to detect gamma-ray clusters in the past x-ray or optical surveys, and our scenario is in accord with apparent no-associations between unidentified EGRET sources and x-ray or optical clusters.

It will be of great significance to perform x-ray or optical observations to search for such loose clusters of galaxies in the regions of high-latitude unidentified EGRET sources. The future gamma-ray projects such as GLAST will also provide a direct test of our scenario. If our scenario is true, a new population of “gamma-ray clusters” will provide us in the future a new probe of dynamically evolving structures in the universe that cannot be traced by x-ray or optical clusters of galaxies.

We would like to thank S. Inoue and T. Naito for useful discussions. TT has partially been supported by the Grant-in-Aid for the Scientific Research Fund (No. 12047233) of the Ministry of Education, Science, and Culture of Japan. TK has been supported in part by the Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists.


format: text

text:
# Statistical geometry in scalar turbulence

\[

## Abstract

A general link between geometry and intermittency in passive scalar turbulence is established. Intermittency is qualitatively traced back to events where tracer particles stay for anomalousy long times in degenerate geometries characterized by strong clustering. The quantitative counterpart is the existence of special functions of particle configurations which are statistically invariant under the flow. These are the statistical integrals of motion controlling the scalar statistics at small scales and responsible for the breaking of scale invariance associated to intermittency.

\]

Scalar fields transported by turbulent flow occur in many physical situations, ranging from the dynamics of the atmosphere and the ocean to chemical engineering (see, e.g., Ref. ). Specific examples are provided by pollutant density, temperature or humidity fields and the concentration of chemical or biological species. The advection-diffusion equation governing the transport of the scalar field $`\theta `$ is :

$$_t\theta (𝒓,t)+\left(𝒗\mathbf{}\right)\theta (𝒓,t)=\kappa \mathrm{\Delta }\theta (𝒓,t),$$
(1)
where $`𝒗(𝒓,t)`$ is the incompressible advecting flow and $`\kappa `$ is the molecular diffusivity. Two broad cases are distinguished : active scalars, where $`𝒗`$ depends on $`\theta `$, e.g. by an explicit relation $`𝒗=𝒗\left(\theta \right)`$, and passive scalars, where the statistics of $`𝒗`$ is independent of $`\theta `$. Here, we shall be concerned with the latter, although we conjecture that the physical mechanisms presented in the following are quite general and relevant also for the active cases. The Fokker-Planck equation (1) is associated to the Lagrangian dynamics of tracer particles whose position $`𝝆(t)`$ obeys $`d𝝆(t)=𝒗(𝝆(t),t)dt+\sqrt{2\kappa }d𝜷(t)`$, where $`𝜷(t)`$ is the isotropic Brownian motion . The equation (1) governs the evolution of the probability density of particles at position $`𝒓`$ and time $`t`$.

Scalar turbulence is typically generated by maintaining a mean scalar gradient $`\theta =𝒈𝒓`$, e.g. by heating/cooling devices in temperature field experiments. The notation $``$ denotes the average with respect to the velocity statistics, which is in principle arbitrary. We shall be interested in flows with correlations having a nontrivial power law behavior in the inertial range of scales $`rL`$, where $`L`$ is the velocity correlation length. Examples are provided by the two and three dimensional Navier-Stokes turbulent flow (see, e.g., Ref. ). A very robust feature of scalar turbulence is its strong intermittency : rare strong events (such as the sharp cliffs observed in Fig. 1) dominate the scalar statistical properties. More quantitatively, intermittency reflects in the anomalous scaling of the correlations. In the inertial range, the scalar structure functions $`S_n(𝒓)=\left(\theta (𝒓,t)\theta (\mathrm{𝟎},t)\right)^n`$ take the form

$$S_n(𝒓)r^{\zeta _n^{dim}}\left(\frac{L}{r}\right)^{\zeta _n^{dim}\zeta _n},$$
(2)
with a nonvanishing value of $`\zeta _n^{dim}\zeta _n`$. Here, $`\zeta _n^{dim}`$ is the value predicted by simple mean field dimensional arguments, e.g. of the Kolmogorov 1941 type. The positive value of the anomalous correction $`\zeta _n^{dim}\zeta _n`$ reflects the breaking of scale-invariance : the scale $`L`$ explicitly appears in the inertial range expressions of scalar correlations, even though $`rL`$. The phenomenon of intermittency is quite generic for scalar turbulence and independent of the specific choice of $`𝒗`$, including for Gaussian velocity fields (see, e.g., Ref. ). Past research on intermittency has mostly concentrated on phenomenological models making vague contacts with the dynamics and based uniquely on scaling exponents and structure functions. Angular dependencies and geometrical informations about the correlations were discarded. It is shown here that the whole structure of multi-point correlations is in fact needed to achieve a real dynamical understanding based on the geometry of the figures identified by tracer particles. Deviations of exponents from naïves dimensional values turn out to be just by-products of the figure geometry nontrivial evolution under Lagrangian dynamics. While the figure size grows according to simple dimensional arguments, the evolution of their shape is a more delicate issue. There are in particular some special functions of the particle positions whose average with respect to the Lagrangian dynamics remains constant in time. That is due to a delicate compensation between the growth due to the figure size and the depletion associated to the figure shape. Those statistically conserved functions are dominating the behavior in the inertial range and controlling the scalar field intermittency.

Specifically, the velocity field $`𝒗`$ considered here is a two-dimensional turbulent flow generated by an inverse energy cascade process . This is a situation of interest both for experiments  and in the atmosphere . The flow realizes the type of turbulence theoretically postulated by A.N. Kolmogorov in 1941 : it is isotropic, it has a constant energy flux (but upscale) and it is scale-invariant with scaling exponent $`1/3`$ . A property of interest to us is that the velocity is not intermittent. All nontrivial scaling properties of the scalar field presented in the following are therefore entirely due to the

advection-diffusion equation (1) and not mere footprints of the velocity field. Details on the integration procedure used for the numerical simulations of (1) with a mean gradient $`𝒈`$ can be found in Ref. . The single-time scalar statistics at the stationary state is defined by the $`n`$-point correlations $`C_n(\underset{¯}{𝒓},t)=\theta (𝒓_1,t)\mathrm{}\theta (𝒓_n,t)`$, where $`\underset{¯}{𝒓}`$ denotes the set $`𝒓_1,\mathrm{},𝒓_n`$. For spatially homogeneous situations, $`C_n`$ is invariant under translations and it depends only on $`2n2`$ degrees of freedom, associated to the separation vectors among the $`n`$ points. In the inertial range, the velocity scale invariance is expected to reflect in scalar correlations $`C_n(\underset{¯}{𝒓})`$ behaving as power laws with respect to global size variables, such as e.g. the gyration radius of the set $`\underset{¯}{𝒓}`$.

Let us consider for simplicity the third order case $`n=3`$. All the following arguments are easily generalized to higher order correlations. The correlation function $`C_3`$ depends on the size, the orientation and the shape of the triangle defined by the three points $`𝒓_1`$, $`𝒓_2`$ and $`𝒓_3`$. The global size variable can be defined as $`R^2=(r_{12}^2+r_{23}^2+r_{31}^2)/3`$, where $`r_{ij}`$ is the distance between the $`i`$-th and the $`j`$-th particle. As shown in Fig. 2, in the inertial range of scales $`C_3`$ depends on $`R`$ as a power law with the exponent $`\zeta _3=1.25`$. The hallmark of intermittency is in the fact that $`\zeta _3`$ is smaller than the dimensional prediction $`5/3`$ (see Ref. ). As for the shape and the orientation of the triangle, we shall use the same Euler parametrization as in Refs. . Defining $`𝝆_1=\left(𝒓_1𝒓_2\right)/\sqrt{2}`$ and

$`𝝆_2=\left(𝒓_1+𝒓_22𝒓_3\right)/\sqrt{6}`$, the shape of the triangle is controlled by the two variables

$$\chi =1/2\mathrm{tan}^1\left[\frac{2𝝆_1𝝆_2}{(\rho _1^2\rho _2^2)}\right];w=2\frac{|𝝆_1\times 𝝆_2|}{R}.$$
(3)
Some of the shapes associated to different values of $`\chi `$ and $`w`$ are shown in Fig. 2. The global orientation of the triangle with respect to the mean gradient direction $`𝒈`$ is defined by the angle $`\phi `$. It is convenient to decompose $`C_3`$ on the orthogonal basis made of $`\mathrm{cos}\left(\mathrm{}\phi \right)`$ and $`\mathrm{sin}\left(\mathrm{}\phi \right)`$. Reversing the coordinates with respect to an axis parallel or orthogonal to $`𝒈`$ statistically leaves the $`\theta `$ field invariant or inverts its sign, respectively. In the projection of $`C_3`$, the angular momentum $`\mathrm{}`$ should therefore be odd and sine functions are absent. Furthermore, the dominant contribution at the small scales is the one having the lowest angular momentum (see Ref.  for the case of Navier-Stokes turbulence). The correlation function $`C_3`$ takes then the form

$$C_3(\underset{¯}{𝒓})=R^{\zeta _3}f(\chi ,w)\mathrm{cos}\phi +\mathrm{},$$
(4)
where the dots stand for subdominant higher-order harmonics of the form $`\mathrm{cos}(2\mathrm{}+1)\phi `$. The invariance under arbitrary permutations of the three vertices of the triangle allows to reduce the phase space to $`\pi /6<\chi <\pi /6`$, $`0<w<1`$ and the function $`f`$ in (4) is antiperiodic in $`\chi `$ with period $`\pi /3`$ . The measured dependence of $`f`$ on the shape coordinates $`\chi `$ and $`w`$ is shown in Fig. 3. The maximum of $`f`$ is realized at $`\chi =0`$, $`w=0`$, where two of the three particles are stuck together. For equilateral triangles ($`w=1`$) or for the “dumbbell” configuration $`\chi =\pi /6`$, $`w=0`$, the symmetries enforce $`f=0`$.

Let us establish the connection with the Lagrangian dynamics. That is done by replacing the Eulerian variables $`\underset{¯}{𝒓}`$ in the argument of $`C_3`$ by their Lagrangian evolutions $`\underset{¯}{𝝆}(t)`$, with $`\underset{¯}{𝝆}(0)=\underset{¯}{𝒓}`$. The resulting object is a

stochastic function whose average with respect to the Lagrangian trajectory statistics is denoted by $`_{}`$. More generally, for a generic function $`\varphi (\underset{¯}{𝒓})`$ of the $`n`$ points $`𝒓_i`$ we can define its Lagrangian average as

$$\varphi (t)_{}=\varphi (\underset{¯}{𝝆})P_n(t,\underset{¯}{𝝆}|0,\underset{¯}{𝒓})𝑑\underset{¯}{𝝆}.$$
(5)
Here, the $`n`$-particle propagator $`P_n(t,\underset{¯}{𝝆}|0,\underset{¯}{𝒓})`$ denotes the probability that, being in $`\underset{¯}{𝒓}`$ at time $`0`$, the $`n`$ particles are in $`\underset{¯}{𝝆}`$ at time $`t`$. In a turbulent Kolmogorov flow, distances typically grow with time as $`|t|^{3/2}`$, whose special instance is the celebrated Richardson law $`r_{12}^2_{}(t)|t|^3`$ for the distance $`r_{12}`$ between two particles. For generic functions $`\varphi `$, homogeneous of positive degree $`\sigma `$ that is $`\varphi (\lambda \underset{¯}{𝒓})=\lambda ^\sigma \varphi (\underset{¯}{𝒓})`$, the Lagrangian average is therefore expected to grow as $`|t|^{3\sigma /2}`$.

Intermittency is dynamically originated by blatant violations of expectations à la Richardson. The anomalous part of correlation functions has in fact a constant Lagrangian average, as clearly demonstrated in Fig. 4 for the third order case. Those statistically preserved functions are the statistical integrals of motion responsible for the breaking of scale invariance associated to intermittency. The other important point is that their constancy is tightly related to the geometry of the figures identified by the Lagrangian particles. Fig. 4 indicates indeed that the size factor $`R^{\zeta _3}`$ grows as $`|t|^{3\zeta _3/2}`$. The Lagrangian average of $`C_3`$ remaining constant, the shape part $`f(\chi ,w)\mathrm{cos}\phi `$ in (4) must compensate for the growth of the figure size. As shown in Fig. 3, the function $`f`$ decreases going from degenerate triangles with two of the vertices close to each other to triangles with aspect ratios of order unity. The geometrical meaning of anomalous scaling laws is then quite clear : the smaller is $`\zeta _3`$, the slower is the compensation needed from the shape factor and the longer degenerate triangle configurations persist. In other words, “stronger intermittency = particles staying longer close to each other”.

Systematic support to the previous physical ideas can be provided in the special case where the advecting velocity in (1) has a short correlation time, the so-called Kraichnan model . The assumption is of course far from realistic, but it leads to the peculiar property that the $`n`$-particle propagators $`P_n`$ in (5) obey closed Fokker-Planck equations  (see also Ref. ). The statistically preserved functions are now identified as zero modes of the $`n`$-particle Fokker-Planck operator . Their anomalous scaling behavior could be calculated in some perturbative limits  or measured numerically . Zero modes enter the $`P_n`$’s via the asymptotic expansion  (see also Ref. ) :

$$P_n(t,\underset{¯}{𝝆}|0,\lambda \underset{¯}{𝒓})=\underset{i,q}{}\lambda ^{\sigma _{i,q}}\varphi _{i,q}(\underset{¯}{𝒓})\psi _{i,q}(t,\underset{¯}{𝝆}),$$
(6)
valid for small $`\lambda `$’s. Zero modes are the $`q=0`$ terms in (6) and they are ordered according to their scaling dimension by the index $`i`$. Higher $`q`$’s identify the so-called slow modes, whose Lagrangian average is growing as an integer power of time, although with an exponent smaller than the dimensional one $`3\sigma _{i,q}/2`$.

A simple example of slow mode for our inverse cascade velocity field is provided for $`n=2`$. The Lagrangian average of $`\left(𝒈𝒓_{12}\right)`$ is preserved as its time derivative is proportional to $`\left(𝒗_1𝒗_2\right)=0`$. The first slow mode associated to it is given by the function $`\left(𝒈𝒓_{12}\right)r_{12}^{2/3}`$. Fig. 5 shows indeed that its Lagrangian average grows as $`|t|`$, much slower at large times than the dimensional law $`|t|^{5/2}`$.

What is the degree of generality of the Lagrangian preservation mechanism for intermittency and the expansion (6) ? The crucial property of the Kraichnan velocity field is its short correlation time (this ensures that the

Lagrangian trajectories are Markovian). For the inverse energy cascade flow considered here this property is lost as the correlation time is finite. The numerical results presented here give therefore a strong indication that the basic mechanisms for scalar intermittency are quite robust and generic. Changing the statistics of the flow affects only quantitative details, such as the numerical value of the exponents or the precise shape of the Lagrangian preserved functions. The expansion (6) is analogous to those encountered in many-body statistical problems  and controls the small-scale behavior of scalar correlations. Indeed, it follows from (1) that the scalar is preserved along the Lagrangian trajectories. Scalar correlations can be then expressed as

$$C_n(\underset{¯}{𝒓},t)=P_n(t,\underset{¯}{𝝆}|0,\underset{¯}{𝒓})\left(𝒈𝝆_1\right)\mathrm{}\left(𝒈𝝆_n\right)𝑑\underset{¯}{𝝆}.$$
(7)
Inserting (6) into this expression it is evident that the behavior at small scales is a superposition of the functions $`\varphi _{i,q}`$. Two general remarks on scalar turbulence follow. First, changing the initial condition $`𝒈𝒓`$ (or the injection mechanism in the forced case) modifies the constants but not the scaling exponents in the correlation functions. This naturally explains the universality properties of scalar turbulence observed in Ref. . Second, all fewer-particle modes $`\varphi _{i,q}`$, i.e. those depending on $`m<n`$ variables, appear in the expression for $`C_n`$. Their Lagrangian averages with respect to the evolution of the $`m`$ particles or the whole set of $`n`$ particles are indeed trivially coinciding. Note however that structure functions satisfy the trivial identity $`S_n(r)=_0^r\mathrm{}_0^r_{r_1}\mathrm{}_{r_n}C_n(\underset{¯}{r})d\underset{¯}{r}`$. All $`m`$-particle modes will therefore drop out from the $`n`$-th order structure function. This clearly illustrates a point previously mentioned : the apparent simplification of the single distance $`r`$ left in structure functions is purely illusory as their anomalous scaling laws are still dynamically associated to $`n`$-particle geometries.

In conclusion, we have identified the integrals of motion responsible for the failure of dimensional analysis and intermittency in scalar turbulence. The anomalous parts of correlation functions are preserved under the Lagrangian dynamics. That preservation is a geometrical effect associated to figures identified by scalar particles persisting in degenerate shapes. Particles tend indeed to remain close to each other for unexpectedly long times, pointing to applications for the reaction rate enhancement in the transport of chemically reactive species.

Acknowledgments We are grateful to U. Frisch and K. Gawȩdzki for illuminating discussions. A. Celani was supported by a “H. Poincaré” CNRS postdoctoral fellowship. The work was partially supported by the EU contract HPRN-CT-2000-00162 “Nonideal Turbulence” and the NSF under Grant No. PHY94-07194.


text:
# Chiral Symmetry, Quark Mass, and Scaling of the Overlap Fermions

## Abstract

The chiral symmetry relation and scaling of the overlap fermions are studied numerically on the quenched lattices at 3 couplings with about the same physical volume. We find that the generalized Gell-Mann-Oakes-Renner relation is satisfied to better than 1% down to the smallest quark mass at $`m_0a=0.006`$. We also obtain the quark mass from the PCAC relation and the pseudoscalar masses. The renormalization group invariant quark mass is shown to be fairly independent of scale. The $`\pi `$ and $`\rho `$ masses at a fixed $`m_\pi /m_\rho `$ ratio indicate small $`O(a^2)`$ corrections. It is found that the critical slowing down sets in abruptly at a very small quark mass close to those of the physical $`u`$ and $`d`$ quarks.

Recent advances in chiral fermion formulation on the lattice hold a great promise in implementing chiral fermions for QCD at finite lattice spacing . There are distinct advantages over the previous formulations. For example, the Wilson fermion breaks chiral symmetry at finite lattice spacing $`a`$ and, therefore, the task of extrapolating the Monte Carlo results to the continuum and chiral limits requires fine tuning and is often difficult. In this case, one could not expect low energy theorems to be reproduced at finite lattice spacing. Similarly, the staggered fermion can only be formulated with four flavors with an $`U(1)`$ subgroup of the flavor non-singlet chiral symmetry. At finite lattice spacing where the numerical calculations are carried out, the flavor symmetry is broken. Furthermore, the anomalous chiral Ward identity does not hold for these fermions unless at the continuum limit. This makes the analysis of anomaly on the lattice rather unclear. There is no unambiguous identification of the fermion zero modes with the topology of the background gauge field. These difficulties have rendered the studies of low energy phenomenology of QCD on the lattice unsettling.

The picture has been altered dramatically with the advent of Neuberger’s overlap fermion  which is derived from the the overlap formalism . All of the above-mentioned impediments can be avoided pending pristine numerical simulation. It is shown to have correct anomaly and exact chiral symmetry on the lattice  and there are no order $`a`$ artifacts . The overlap fermion has a compact form in four dimension and is easily employed to derive low energy theorems, chiral symmetry relations, and anomaly at finite lattice spacing. In this letter, we shall study the overlap fermion numerically. We test chiral symmetry via the Gell-Mann-Oakes-Renner relation at finite lattice spacing and check the scaling behavior of the $`\pi `$ and $`\rho `$ masses. The preliminary results were reported earlier . We also obtain the quark mass from the Chiral Ward identity and verify that it is free from additive renormalization and the renormalization group invariant quark mass is indeed independent of scale.

Neuberger’s Dirac operator has the following form for the massive case

$$D(m_0)=1+\frac{m_0a}{2}+(1\frac{m_0a}{2})\gamma _5ϵ(H),$$
(1)
where $`ϵ(H)=H/\sqrt{H^2}`$ is the matrix sign function of H which we take to be the Hermitian Wilson-Dirac operator, i.e. $`H=\gamma _5D_w`$. Here $`D_w`$ is the usual Wilson fermion operator, except with a negative mass parameter which corresponds to $`\kappa _c<\kappa <0.25`$. We take $`\kappa =0.19`$ in our calculation. The massless operator $`D(0)`$ is shown  to satisfy the Ginsparg-Wilson relation : $`\{\gamma _5,D(0)\}=D(0)\gamma _5D(0)`$. The bare mass parameter $`m_0`$ is proportional to the quark mass without additive constant which we shall verify.

There are several numerical approaches to approximate the sign function $`ϵ(z)`$ . We adopt the optimal rational approximation  with a ratio of polynomials of degree 12 in the Remez algorithm. We find the error to the approximation of $`ϵ(z)`$ to be within $`10^5`$ in the range \[0.02, 2\] of the argument $`z`$. In the range \[0.0005, 0.02\], the error can be as large as 1%. To improve the accuracy of $`ϵ(H)`$ and hence the chiral symmetry property, the smallest 10 to 20 eigenvalues of $`H^2`$ with eigenvalues of $`|H|`$ less than 0.04 are projected out for exact evaluation of the sign function from these eigenstates . This has the added benefit of reducing the number of iterations for the multi-mass conjugate gradient inversion of $`H^2+c_i`$ in the inner loop by a factor of 3.5 or so. We checked the unitarity of the matrix $`V=\gamma _5ϵ(H)`$. For $`Vx=b`$, we find $`|x^{}xb^{}b|10^9`$. Since $`V`$ is unitary we exploit the identity $`(1+V^{})(1+V)=2+V+V^{}`$ in order to use the conjugate gradient algorithm on the hermitian matrix $`V+V^{}`$ instead of $`V^{}V`$ which has a higher condition number. Furthermore, since $`[V+V^{},\gamma _5]=0`$, one can use a chiral source, i.e. $`\gamma _5b=\pm b`$ to save one matrix multiplication  per iteration. In the present study, we consider three lattices: $`6^3\times 12`$ at $`\beta =5.7`$, $`8^3\times 16`$ at $`\beta =5.85`$, and $`10^3\times 20`$ at $`\beta =6.0`$, which have about the same physical volume. With residuals at $`10^7`$, the inner loop takes typically $`200250`$ iterations and is almost independent of the lattice volume. The outer loop takes $`40110`$ iterations from the small to large lattice volumes. Even for topological sectors with charge $`Q0`$, we find the critical slowing down is much milder than that of the Wilson fermion and we found no exceptional configurations. The critical slowing down sets in quite abruptly when $`m_0a<0.006`$ for the sector with topology. This is already very close to the physical $`u`$ and $`d`$ masses.

The lattice chiral symmetry is reflected in the generalized Gell-Mann-Oakes-Renner (GOR) relation

$$m_0ad^4x\pi (x)\pi (0)=2\overline{\psi }\psi .$$
(2)
Since it is satisfied for each quark mass, lattice volume and spacing, configuration by configuration, and for each source , it serves as an economic test of chiral symmetry. Here $`\pi (x)=\overline{\psi }\gamma _5(\tau /2)\psi `$ is the pion interpolation field and the quark propagator is the external one  with $`D_c^1=(1m_0a/2)^1[D^1(m_0)1/2]`$. Alternatively, one can use the bilinears $`(1m_0a/2)^1\overline{\psi }\mathrm{\Gamma }(1D/2)\psi `$ for the operators and $`D^1`$ as the propagator. We utilize this relation in Eq. (2) as a check of the numerical implementation of the Neuberger operator. We find that for the lattices we consider the GOR relation is satisfied very well (to within $`10^3`$) all the way down to the smallest mass $`m_0a=2\times 10^4`$ for the $`Q=0`$ sector. For the $`Q0`$ sector, the presence of zero modes demands higher precision for the approximation of $`ϵ(H)`$.

For example, we show in Fig. 1 the ratio of the right to left side of Eq. (2) for a typical configuration with topology on the $`6^3\times 12`$ lattice at $`\beta =5.7`$ as a function of the quark mass $`m_0a`$. When only 10 smallest eigenmodes of $`H^2`$ are projected, we see that the ratio deviates from one for small quark masses. For the smallest mass $`m_0a=0.006`$, it may deviate as much as 16%. The situation is considerably improved when 20 smallest eigenmodes are included where the deviation is reduced to 1%. Our result with the overlap is appreciably better than the domain-wall fermion case when the size of the fifth dimension is limited to $`L_s=10`$ to 48 . In the latter case, the ratio deviates from unity by $`55\%`$ for $`L_s=10`$ and $`15\%`$ for $`L_s=48`$ at the quark mass $`m_fa=0.02`$ which is about 3 times larger than the mass of $`m_0a=0.006`$ in our study where the deviation is at the 1% level.

The average of $`u`$ and $`d`$ quark masses in the $`\overline{MS}`$ scheme at the renormalization scale $`\mu `$ is obtained from the axial Ward identity $`Z_A_\mu A_\mu =2Z_S^1m_0Z_PP`$ via the ratio of the matrix elements

$$m_q^{\overline{MS}}(\mu )=Z_S^1(\mu )m_0=\underset{t\mathrm{}}{lim}\frac{Z_A(\mu )}{Z_P(\mu )}\frac{_\stackrel{}{x}0|_tA_4(x)|\pi (0)}{2_\stackrel{}{x}0|P(x)|\pi (0)}$$
(3)
where $`A_\mu =\overline{\psi }i\gamma _\mu \gamma _5(\tau /2)\psi `$ and $`P=\overline{\psi }i\gamma _5(\tau /2)\psi `$. The resumed cactus diagram of one-loop calculation of $`Z_A`$ and $`Z_P`$  are used in Eq. (3). It turns out that the ratio $`Z_A/Z_P`$ is very close to unity for $`\mu =1/a`$. The renormalized quark mass $`m_q^{\overline{MS}}a`$ defined in Eq. (3) is plotted in Fig. 2 against the bare mass $`m_0a`$ for the three lattices.

We first observe that the renormalized quark mass does not have an additive part due to the lattice chiral symmetry. The linear fit including the smallest 5, 7, and 8 quark masses for $`\beta =5.7,5.85`$, and 6.0 respectively shows that the intercepts are of the order of $`10^3`$ and are consistent with zero. This is a distinct advantage over the Wilson fermion whose quark mass is subject to additive renormalization which depends on the gauge configuration. From the slope we can determine the non-perturbative $`Z_S^1`$ which are 0.95(5), 0.89(5) and 0.80(1) for $`\beta =5.7,5.85`$, and 6.0. These are within 14% from the perturbative calculation which are 0.818, 0.824, and 0.829 respectively . The renormalization group invariant quark mass is defined as the integration constant of the evolution equation such that

$$m_q^{RGI}=\mathrm{\Delta }Z_S(\mu )m_q^{\overline{MS}}(\mu )=\mathrm{\Delta }Z_S(\mu )Z_S^1(\mu )m_0.$$
(4)
Using the four-loop calculation of $`\mathrm{\Delta }Z_S(\mu )`$ in the $`\overline{MS}`$ scheme , we obtain the product of $`\mathrm{\Delta }Z_S(\mu )Z_S^1(\mu )`$ which should be scale invariant. The results from the 3 lattices give 1.12(1), 1.17(6), and 1.14(5) for $`\beta =5.7,5.85`$, and 6.0. They are indeed independent of scale within errors. From the pion mass fit to be discussed later, we determine $`m_0=6.8(5)`$ MeV from the physical pion mass on the $`\beta =6.0`$ lattice. This gives $`m_q^{RGI}=7.6(6)`$ MeV for the average of $`u`$ and $`d`$ quark masses, which in turn gives $`m_q^{\overline{MS}}(\mu =2\mathrm{G}\mathrm{e}\mathrm{V})=5.5(5)`$ MeV. The lattice scale is set by the $`r_0`$ as derived from the static quark potential . We should point out that this quark mass should not be taken literally, since the volume is quite small. The primary purpose of the present study is to verify that the quark mass extracted from the overlap fermion action does not have additive renormalization and that $`m_q^{RGI}`$ is indeed scale invariant.

Next we look at hadron masses. The results for the $`\pi ,\rho `$ and nucleon masses on the $`10^3\times 20`$ lattice at $`\beta =6.0`$ are given in Fig. 3. They are obtained by a single exponential fit with covariance matrix.

We see clearly that the nucleon mass suffers from the finite volume effect when $`m_0`$ is smaller than 0.4 GeV. Although the $`\rho `$ and $`\pi `$ masses appear not affected as much, there is the worry that the finite volume effect is present when one considers the pion mass behavior as a function of the quark mass $`m_0`$. In the finite size scaling region where the lattice size is much smaller than the pion Compton wavelength ($`L<<1/m_\pi `$), the eigenvalue distribution and the chiral condensate can be described by the universal function derived in the chiral random matrix theory . While the chiral condensate of the overlap fermion in this region has been studied , we shall concentrate on the region $`L>1/m_\pi `$ where the chiral perturbation analysis is expected to apply. Plotted in Fig. 4 are the pseudoscalar meson mass squared ($`m_P^2a^2`$) as a function of $`m_0a`$ for the three lattices for those pseudoscalar mesons whose Compton wavelengths are less than the respective lattice size.

We fit them with a form suggested by the quenched chiral perturbation theory

$$m_P^2a^2=2Am_0a^2\{1\delta \mathrm{ln}(2Am_0/\mathrm{\Lambda }_\chi ^2)\}+4Bm_0^2.$$
(5)
We find that for $`\mathrm{\Lambda }_\chi `$ in between 0.6 GeV and 1.4 GeV, the chiral log $`\delta `$ for the $`\beta =6.0`$ case is in the range of 0.20(2) to 042(14) which is slightly larger than what is obtained from the Wilson fermion . The fit is quite stable with small $`\chi ^2`$ ($`\chi ^2/DF<0.1`$) and is fairly insensitive to the range of the fitted quark masses shown in Fig. 4. On the other hand, the fit for $`\beta =5.85`$ and 5.7 lattices shows that $`\delta `$ is consistent with zero (the typical value ranges from -0.04(21) to -0.34(48)) and the result is more sensitive to the range of quark masses that are fitted. We are not able to draw any definite conclusion from these inconsistent results, except to point out that perhaps this is due to the finite volume effect. After all, the physical volumes of the lattices we use are quite small ($`L1fm`$). Another possibility is that the smaller quark mass cases may already fall into the finite size scaling region where the pion mass may have a different behavior than prescribed in Eq. (5). This will be studied elsewhere.

Finally, we examine the scaling of $`\rho `$ and $`\pi `$ masses. Since $`m_\pi ^2a^2`$ and $`m_\rho a`$ are fairly linear in $`m_0a`$, we fit them with $`m_\rho a=A+Bm_0a`$ and $`m_\pi ^2a^2=Cm_0a+Dm_0^2a^2`$ for simplicity. The fits are decent with $`\chi ^2/DF<1`$ for the three lattices and the full range of quark masses in Fig. 4. From the fits of $`m_\rho a`$ and $`m_\pi ^2a^2`$, we determine $`m_0`$ for which the ratio of $`m_\pi /m_\rho =0.4,0.5,0.6`$ and plot in Fig. 5 the corresponding $`m_\rho `$ and $`m_\pi `$ in units of $`\sqrt{\sigma }`$ as a function of $`\sigma a^2`$, where $`\sigma `$ is the string tension. The errors on the vector and pseudoscalar masses are determined from interpolating the data from the neighboring quark masses. It is known that the overlap operator does not have $`O(a)`$ artifacts . It appears from Fig. 5 that even the $`O(a^2)`$ errors are small.

To conclude, we find that when the matrix sign function $`ϵ(H)`$ in the overlap fermion is well approximated, the promised chiral symmetry at finite lattice spacing and the scaling of the renormalization group invariant quark mass and hadron masses are manifested in the present numerical calculation. One drawback of the overlap fermion is the large numerical overhead in the present algorithm. But, the stake of being able to implement chiral symmetry at finite lattice spacing is high. Furthermore, the unexpected feature of being able to push the critical slowing down to close to the physical $`u`$ and $`d`$ quark masses has the advantage of being able to study the correct chiral behavior. The nice scaling result is encouraging for controlling the continuum extrapolation and may afford the possibility of working at relatively large lattice spacings. For the moment, the study is limited to small volumes. As long as one can extend it to large volumes, it appears that one will be at a stage of putting all the systematic errors of at least the quenched approximation under control.

This work is partially supported by DOE Grants DE-FG05-84ER40154 and DE-FG02-95ER40907. We thank R. Edwards for sharing his experience in implementing the sign function solver. We also thank H. Neuberger and T.W. Chiu for many stimulating discussions, R. Narayanan for communication, and H. Panagopoulos for calculating the renormalization constants for us.