text:
# II. Territory covered by 𝑁 random walkers on stochastic fractals. The percolation aggregate

## I INTRODUCTION

In the preceding paper  we studied the problem of evaluating the average number $`S_N(t)`$ of distinct sites (or territory covered) by $`N`$ independent random walkers that, all starting from the same site, diffuse on a deterministic fractal during time $`t`$. An answer to this problem was given in Ref.  in terms of an asymptotic series completely analogous to that given in Refs.  for Euclidean substrates. In this present paper we turn our attention to the case of stochastic (or disordered) fractal media and again propose the same kind of asymptotic solution, although, in contrast with all the previous cases, we now have to translate the procedure to the topological or chemical language.

Stochastic media are not constructed by the iteration of an invariable rule, such as that used in deterministic fractals, but are rather the result of a random process. Consequently, the concept of generator and rigorous self-similarity is absent and their fractal nature is recognized by the scaling of statistical quantities. Many natural objects share this statistical-fractal structure  so that stochastic models seem to be more suitable to represent diffusion in real media. Of particular interest is the percolation model that has been used to characterize many disordered systems . This model is constructed by filling a regular lattice with “occupied” sites with a certain probability. Nearest neighbor occupied sites are supposed to be connected and form a series of clusters. At a certain critical concentration $`p_c`$ an infinite cluster appears, which is called the incipient percolation aggregate or percolation cluster at criticality.

The quantity we are interested in, $`S_N(t)`$, is, for disordered media, the result of a double average: an average over the walks that the $`N`$ random walkers can perform over a given lattice, followed by an average over many (ideally, all) realizations of the random lattice. This situation leads to certain subtleties, absent in deterministic fractals, that require special treatment. In particular, $`S_N(t)`$ can be expressed by

$$S_N(t)=\left\{1\left[\mathrm{\Gamma }_t(𝐫)\right]^N\right\},$$
(1)
where the sum is over all the sites of each fractal lattice, $`\mathrm{\Gamma }_t(𝐫)`$ is the survival probability, i.e., the probability that site $`𝐫`$ has not been visited by time $`t`$ by a single random walker starting from the origin, $`\left\{1\left[\mathrm{\Gamma }_t(𝐫)\right]^N\right\}`$ represents the mean territory explored by the $`N`$ random walkers on a given lattice (the first average) , and $`[\mathrm{}]`$ indicates that the average (the second average) of $`[\mathrm{}]`$ has to be performed over all possible stochastic lattices compatible with the random generation rules. Equation (1) can be rewritten as

$$S_N(t)=\underset{m=0}{\overset{\mathrm{}}{}}\underset{i=1}{\overset{n(m)}{}}\left\{1\left[\mathrm{\Gamma }_t(𝐫_{m,i})\right]^N\right\},$$
(2)
where $`𝐫_{m,i}`$ stands for the $`i`$-th site out of $`n(m)`$ that are separated from the origin by a Euclidean distance between $`m\mathrm{\Delta }rr_m`$ and $`(m+1)\mathrm{\Delta }r`$ with $`\mathrm{\Delta }r`$ small (say, of the order of the lattice spacing). If $`\mathrm{\Gamma }_t(𝐫_{m,i})`$ is almost independent of $`i`$ and the lattice realization, i.e., if the fluctuations in the probability density $`\mathrm{\Gamma }_t(𝐫_{m,i})`$ follow a narrow distribution, then one could approximate $`\mathrm{\Gamma }_t(𝐫_{m,i})\mathrm{\Gamma }_t(𝐫_{m,i})\mathrm{\Gamma }_t(r_m)`$, and therefore estimate $`S_N(t)`$ by

$$S_N(t)=\underset{m=0}{\overset{\mathrm{}}{}}\left\{1\left[\mathrm{\Gamma }_t(r_m)\right]^N\right\}n(m),$$
(3)
where $`n(m)`$ is the average number of fractal sites separated from the origin by a distance bracketed by $`r_m`$ and $`r_m+\mathrm{\Delta }r`$. This is essentially the starting relationship used (implicitly) by Havlin et al.  to find that, for large $`N`$,

$$S_N(t)t^{d_s/2}(\mathrm{ln}N)^{d_f/u}$$
(4)
in the non-trivial time regime (or regime II) . Here $`u=d_w/(d_w1)`$, $`d_w`$ is the anomalous diffusion exponent, $`d_f`$ is the fractal dimension of the substrate, and $`d_s=2d_f/d_w`$ is the spectral dimension. However, the hypothesis leading to Eq. (3) is in general false, as we will explicitly show in Sec. II by means of numerical simulations of $`\mathrm{\Gamma }_t(𝐫_{m,i})`$ for the two-dimensional percolation cluster at criticality. Indeed, it is known  that the fluctuations of the probability density $`P(r,t)`$ of random walks (also called the propagator or Green’s function), which is a statistical quantity closely related to the survival probability, exhibits a broad logarithmic distribution for some random fractals such as percolation clusters and self-avoiding walks. Bunde et al.  have found that the quantity $`P(r,t)^q`$ exhibits multifractal scaling, $`P(r,t)^qP(r,t)^{\tau (q)}`$, where $`\tau (q)q^\gamma `$ and $`\gamma =(d_w^{\mathrm{}}1)/(d_w1)`$ and $`d_w^{\mathrm{}}`$ is the chemical random walk dimension. This behavior is a consequence of the large fluctuations of $`P(r,t)`$ for fixed $`r`$ and $`t`$ from a given aggregate to another. Nevertheless, these authors have also shown that the distribution of the propagator in the chemical $`\mathrm{}`$ space, $`P(\mathrm{},t)`$, is narrow and, consequently, $`P(\mathrm{},t)^qP(\mathrm{},t)^q`$. The chemical distance $`\mathrm{}`$, the length of the shortest path between two sites along lattice bonds, is a more natural measure than the Euclidean distance in disordered systems. It is, for example, the distance used in the calculation of optimum paths in cities.

Let $`\mathrm{}_{m,i}`$ label the $`i`$-th site out of those $`n(m)`$ that are placed at a chemical distance $`\mathrm{}`$ from a given origin with $`\mathrm{}_m\mathrm{}<\mathrm{}_{m+1}`$, $`\mathrm{}_m=m\mathrm{\Delta }\mathrm{}`$ and $`\mathrm{\Delta }\mathrm{}`$ small (say, of the order of the lattice spacing), and let $`\mathrm{\Gamma }_t(\mathrm{}_{m,i})`$ be the survival probability in the chemical space defined as the probability that site $`\mathrm{}_{m,i}`$ has not been visited by time $`t`$ by a single random walker starting from the origin. Then we can rewrite Eq. (2) in the chemical $`\mathrm{}`$ space as

$$S_N(t)=\underset{m=0}{\overset{\mathrm{}}{}}\underset{i=1}{\overset{n(m)}{}}\left\{1\left[\mathrm{\Gamma }_t(\mathrm{}_{m,i})\right]^N\right\}.$$
(5)
One may expect that the distribution of $`\mathrm{\Gamma }_t(\mathrm{}_{m,i})`$ for fixed $`\mathrm{}_m`$ and $`t`$ is as narrow as the distribution of the propagator in the chemical space. (We support this conjecture in Sec. III by means of numerical simulations of $`\mathrm{\Gamma }_t(\mathrm{}_{m,i})`$ in two-dimensional percolation clusters at criticality.) In this case

$$\mathrm{\Gamma }_t(\mathrm{}_{m,i})\mathrm{\Gamma }_t(\mathrm{}_{m,i})\mathrm{\Gamma }_t(\mathrm{}_m)$$
(6)
for all possible lattice realizations so that $`\left[\mathrm{\Gamma }_t(\mathrm{}_{m,i})\right]^N\mathrm{\Gamma }_t(\mathrm{}_{m,i})^N`$ and, therefore, $`S_N(t)`$ can be approximated by

$$S_N(t)=\underset{m=0}{\overset{\mathrm{}}{}}\left\{1\left[\mathrm{\Gamma }_t(\mathrm{}_m)\right]^N\right\}n(m),$$
(7)
where $`n(m)`$ is the average number of fractal sites separated from the origin by a chemical distance with value between $`\mathrm{}_m`$ and $`\mathrm{}_m+\mathrm{\Delta }\mathrm{}`$. From this formula and following the procedure outlined in the preceding paper , in Sec. II we will arrive at an expression for $`S_N(t)`$ for the non-trivial time regime whose leading asymptotic behavior coincides, apart from the value of the prefactor, with the recent proposal, based on a scaling approach, of Dräger and Klafter :

$$S_N(t)t^{d_s/2}(\mathrm{ln}N)^{d_{\mathrm{}}/v}$$
(8)
with $`v=d_w^{\mathrm{}}/(d_w^{\mathrm{}}1)`$. Equation (8) differs from the relationship proposed by Havlin et al. , Eq. (4), for those cases such as that considered in this present paper where $`d_{\text{min}}1`$. Both Havlin et al. and Dräger and Klafter supported their conjectures by means of data collapsing plots of computer simulation results obtained for two- and three-dimensional percolation aggregates, respectively. We here want also draw attention to the risk involved in this method of analysis when the influence of the corrective terms is not properly considered since these terms have a large influence  on the final value of $`S_N(t)`$.

The paper is organized as follows. In Sec. II we present the asymptotic evaluation of $`S_N(t)`$ for stochastic fractal lattices as a translation to the chemical space of the procedure implemented in the Euclidean and deterministic fractal cases . A less rigorous but fairly simple and instructive method for obtaining the main asymptotic term and estimating the corrective terms of $`S_N(t)`$ for large $`N`$ is also presented. In Sec. III we report simulation results for the survival probability of a random walker on a two-dimensional incipient percolation aggregate when a trap is placed at a site at a fixed chemical distance $`\mathrm{}`$ or Euclidean distance $`r`$. We find that the distribution is narrow \[broad\] if the traps are located at a fixed chemical \[Euclidean\] distance. The parameters governing (i) the asymptotic behavior of $`\mathrm{\Gamma }_t(\mathrm{})`$ ($`c`$, $`v`$, $`\mu `$ and $`A`$), (ii) how the fractal volume grows ($`V_0^{\mathrm{}}`$ and $`d_{\mathrm{}}`$), and (iii) how fast a single walker diffuses ($`D_{\mathrm{}}`$ and $`d_w^{\mathrm{}}`$) are estimated in this section. In Sec. IV we compare the zeroth- and first-order asymptotic expansion for $`S_N(t)`$ with simulation results obtained for the two-dimensional incipient percolation aggregate. We also criticize the reliability of typical collapsing plots for the determination of the dominant trend of $`S_N(t)`$. We conclude with some remarks in Sec. V.

## II TERRITORY COVERED BY $`N`$ RANDOM WALKERS ON A STOCHASTIC FRACTAL SUBSTRATE

In this section we will translate the results of the previous paper  to chemical language. Reasons for this procedure have already been given in Sec. I. We start by replacing Eq. (7) by its continuum approximation

$`S_N(t)`$ $``$ $`{\displaystyle _0^{\mathrm{}}}\left\{1\left[\mathrm{\Gamma }_t(\mathrm{})\right]^N\right\}d_{\mathrm{}}V_0^{\mathrm{}}\mathrm{}^{d_{\mathrm{}}1}𝑑\mathrm{}`$ (9)
$``$ $`V_0^{\mathrm{}}d_{\mathrm{}}(2D_{\mathrm{}})^{d_{\mathrm{}}/2}t^{d_{\mathrm{}}/d_w}{\displaystyle _0^{\mathrm{}}}\left\{1\left[\mathrm{\Gamma }_t(\xi )\right]^N\right\}\xi ^{d_{\mathrm{}}1}𝑑\xi ,`$ (10)
where $`dV(\mathrm{})=V_0^{\mathrm{}}d_{\mathrm{}}\mathrm{}^{d_{\mathrm{}}1}d\mathrm{}`$ is the average number of fractal sites placed at a chemical distance between $`\mathrm{}`$ and $`\mathrm{}+d\mathrm{}`$, and $`\xi \mathrm{}/(\sqrt{2D_{\mathrm{}}}t^{1/d_w^{\mathrm{}}})`$. Here $`D_{\mathrm{}}`$ is the diffusion constant defined by the Einstein relation,

$$L^2=2D_{\mathrm{}}t^{2/d_w^{\mathrm{}}},$$
(11)
where $`L^2\mathrm{}^2`$ is the mean-square chemical distance traveled by a single random walker by time $`t`$ ($`t`$ large). Next, we assume that the asymptotic dependence of $`\mathrm{\Gamma }_t(\mathrm{})`$ for $`\xi 1`$ is given by

$$\mathrm{\Gamma }_t(\mathrm{})1A\xi ^{\mu v}e^{c\xi ^v}\left(1+\underset{n=1}{\overset{\mathrm{}}{}}h_n\xi ^{nv}\right),$$
(12)
with $`v=d_w^{\mathrm{}}/(d_w^{\mathrm{}}1)`$. This functional form holds (at least in its first terms) on Euclidean lattices  and agrees with the expression conjectured in Refs.  for fractal substrates with $`d_{\text{min}}=1`$ such as the Sierpinski gaskets. The dominant asymptotic behavior of the propagator in chemical space , $`P(l,t)\mathrm{exp}\left(c\xi ^v\right)`$, also coincides with the assumed dominant exponential decay of the mortality function $`1\mathrm{\Gamma }_t(\mathrm{})`$ in Eq. (12). It is known that both the propagator and the mortality function share the same asymptotic behavior for Euclidean lattices and for the Sierpinski lattice  and we can expect that this coincidence also is the case for stochastic fractals (we will check this supposition in Sec. III). The rest of the analysis is identical with that carried out for Euclidean and deterministic fractal lattices  save for the change of the Euclidean parameters to their chemical analogues. The result for $`S_N(t)`$ is, consequently,

$$S_N(t)\widehat{S}_N(t)\left(1\frac{d_{\mathrm{}}}{v}\underset{n=1}{\overset{\mathrm{}}{}}\mathrm{ln}^nN\underset{m=0}{\overset{n}{}}s_m^{(n)}\mathrm{ln}^m\mathrm{ln}N\right)$$
(13)
with

$$\widehat{S}_N(t)=V_0^{\mathrm{}}(2D_{\mathrm{}})^{d_{\mathrm{}}/2}t^{d_{\mathrm{}}/d_w^{\mathrm{}}}\left(\frac{\mathrm{ln}N}{c}\right)^{d_{\mathrm{}}/v}$$
(14)
and

$`s_0^{(1)}`$ $`=`$ $`\omega `$ (15)
$`s_1^{(1)}`$ $`=`$ $`\mu `$ (16)
$`s_0^{(2)}`$ $`=`$ $`(\beta 1)\left({\displaystyle \frac{\pi ^2}{12}}+{\displaystyle \frac{\omega ^2}{2}}\right)(ch_1\mu \omega )`$ (17)
$`s_1^{(2)}`$ $`=`$ $`\mu ^2+(\beta 1)\mu \omega `$ (18)
$`s_2^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\beta 1)\mu ^2.`$ (19)
Here $`\omega =\gamma +\mathrm{ln}A+\mu \mathrm{ln}c`$, $`\gamma 0.577215`$ is the Euler constant, and $`\beta =d_{\mathrm{}}/v`$. The dependence on $`t`$ and $`N`$ of the main term of $`S_N(t)`$ as given by Eq. (14), i.e., $`S_N(t)t^{d_s/2}(\mathrm{ln}N)^{d_{\mathrm{}}/v}`$, coincides with the prediction of Ref. .

### A simpler way to estimate the territory covered

We finish this section by showing how to find the full main term of Eq. (13) and even predict the form of the corrective terms by only resorting to extremely simple arguments. The ideas of the following reasoning were already used in Ref. . The crucial point in our argument is that, for a fixed time $`t`$, $`1[\mathrm{\Gamma }_t(\mathrm{})]^N`$ approaches a unit step function $`\mathrm{\Theta }(\mathrm{}\mathrm{}_\times )`$ when $`N\mathrm{}`$, $`\mathrm{}_\times `$ being the step’s width (see Fig. 1). The reason for this behavior is clear: For large $`N`$, $`[\mathrm{\Gamma }_t(\mathrm{})]^N`$ is only non-neglible when $`\mathrm{\Gamma }_t(\mathrm{})`$ is very close to 1. Obviously this occurs when the root-mean square chemical distance $`L(t)`$ traveled by the single random walker by time $`t`$ is small compared with $`\mathrm{}`$, i.e., when $`\xi =\mathrm{}/L(t)`$ is large. This in turn implies that in the evaluation of $`S_N(t)`$ only the behavior of $`\mathrm{\Gamma }_t(\mathrm{})`$ for large $`\xi `$ is relevant. Then, as $`1[\mathrm{\Gamma }_t(\mathrm{})]^N`$ approaches a step function of width $`\mathrm{}_\times `$, the integration of Eq. (9) yields

$$S_N(t)V_0^{\mathrm{}}\mathrm{}_\times ^d_{\mathrm{}},$$
(20)
i.e., the territory covered is just the volume of a chemical hypersphere of radius $`\mathrm{}_\times `$. Defining the width $`\mathrm{}_\times `$ of the step function as the distance at which $`1[\mathrm{\Gamma }_t(\mathrm{})]^N`$ takes the intermediate value $`1/2`$ (any other value between 0 and 1 would also be valid as $`\mathrm{}_\times `$ is not very sensitive to this value when $`N1`$), and assuming that $`1\mathrm{\Gamma }_t(\mathrm{})A\xi ^{\mu v}\mathrm{exp}(c\xi ^v)`$ for large $`\xi `$, we deduce that $`1/2NA\xi _\times ^{\mu v}\mathrm{exp}(c\xi _\times ^v)`$, with $`\xi _\times \mathrm{}_\times /L`$ so that

$$c\xi _\times ^v\mathrm{ln}N\mu v\mathrm{ln}\xi _\times +\mathrm{ln}2A.$$
(21)
The term $`\mathrm{ln}N`$ is dominant on the right-hand side of Eq. (21) for large $`N`$, so that a first-order solution of this equation is

$$c\xi _\times ^v\mathrm{ln}N,$$
(22)
i.e., $`\mathrm{}_\times ^vL^v\mathrm{ln}(N)/c`$. Hence Eq. (20) yields

$$S_N(t)V_0^{\mathrm{}}L^d_{\mathrm{}}\left(\frac{\mathrm{ln}N}{c}\right)^{d_{\mathrm{}}/v}$$
(23)
which is in full agreement with the main term of Eq. (13) when the Einstein relation, Eq. (11), is considered. Inserting the above first-order solution for $`\mathrm{}_\times `$ into the right-hand side of Eq. (21), we get the improved solution $`c\xi _\times ^v\mathrm{ln}N\mu \mathrm{ln}\mathrm{ln}N+\mathrm{ln}Ac^\mu +\mathrm{ln}2`$, so that Eq. (20) becomes

$`S_N(t)`$ $``$ $`V_0^{\mathrm{}}L^d_{\mathrm{}}\left({\displaystyle \frac{\mathrm{ln}N}{c}}\right)^{d_{\mathrm{}}/v}\left[1+{\displaystyle \frac{d_{\mathrm{}}}{v}}{\displaystyle \frac{\mu \mathrm{ln}\mathrm{ln}N+\mathrm{ln}Ac^\mu +\mathrm{ln}2}{\mathrm{ln}N}}\right].`$ (24)
This expression is strikingly close to the first-order approximation of Eq. (13), the only difference being that the term $`\mathrm{ln}2=0.693\mathrm{}`$ in Eq. (24) plays the role of the Euler constant $`\gamma 0.577215`$ in Eq. (13). Finally, notice that this simple method is not limited to disordered media but that it is also valid for estimating $`S_N(t)`$ for the non-disordered substrates (Euclidean and deterministic fractal media) considered in Refs. .

## III SURVIVAL PROBABILITY, FRACTAL VOLUME AND DIFFUSION IN PERCOLATION AGGREGATES

We have carried out simulations for the number of distinct sites visited by $`N`$ independent random walkers on a typical stochastic fractal: the percolation aggregate embedded in two dimensions. In our simulations every random walker makes a jump from a site to one of its nearest neighbors placed at one unit distance in each unit time. The incipient percolation aggregates were constructed by the standard Leath method  on a square lattice with side $`400`$. In the Leath method, a seed is placed on the site in the center of this box and the cluster is generated by epidemic spreading to their nearest neighbors with an infecting probability $`p_c=0.5927460`$  corresponding to site percolation in the square lattice. At every step of the generation process a new chemical shell is added to the previous shell by the occupation with probability $`p_c`$ of its empty nearest neighbors. The process continues until we reach a shell whose sites do not infect any of their empty neighbors; in this case the aggregate so generated is rejected. If the cluster grown in this way spans the box from a side to the other, the cluster is accepted as a good representation of a portion of an infinite percolation aggregate. Our simulations were carried out over 2000 aggregates generated in this way.

In order to compare the simulation results for $`S_N(t)`$ with the predictions of our theoretical approach, Eq. (13), we must check that the survival probability or, equivalently, the mortality function, $`h(\mathrm{},t)=1\mathrm{\Gamma }_t(\mathrm{})`$, really behaves in the form conjectured in Eq. (12). Moreover, we must confirm first that, for a given chemical distance $`\mathrm{}`$, the distribution of $`h(\mathrm{},t)`$ over different realizations of the incipient percolation cluster is narrow because our theoretical analysis \[cf. Eq. (7)\] relies on this assumption \[cf. Eq. (6)\]. The numerical evaluation of this quantity as well as the propagator $`P(𝐢,t)`$ (i.e., the probability of finding a single random walker at site $`𝐢`$ at time $`t`$) is performed by the Chapman-Kolmogorov method (also called the exact enumeration method ). Initially the propagator takes the value $`1`$ at the origin site, $`P(\mathrm{𝟎},0)=1`$, and $`0`$ at any other site. This density evolves at every time step by updating its value at every site in the form described by the master equation

$$P(𝐢,t+1)=\frac{1}{b}\underset{\text{neighbor}=1}{\overset{M}{}}P(𝐣,t)+\left(1\frac{M}{b}\right)P(𝐢,t),$$
(25)
where $`b`$ is the maximum coordination number ($`b=4`$ for the two-dimensional percolation aggregate) and $`M`$ is the number of neighbors of site $`𝐢`$ that belong to the aggregate. This evolution equation is valid for the so-called blind “ants”  because the random walker “attempts” a jump at time $`t+1`$ to a possible nearest neighbor (selected randomly) of the site it occupied at time $`t`$. If the selected site does not belong to the cluster, the random walker stays at the same site and no jump takes place. The “blindness” of the random walkers is taken into account by the second term of the right-hand side of Eq. (25). The trap is simulated by a special site belonging to the cluster that absorbs all the probability density that enters it without giving back any probability to its neighbors. In the simulation of the mortality function, we located a trap at a chemical distance $`\mathrm{}=30`$ in each of the 2000 percolating clusters. We repeated the experiment for traps located at a fixed Euclidean distance, $`r=80`$. The resulting histogram for $`t=1000`$ is shown in Fig. 2 (to be compared with the histogram of the propagator shown in Fig. $`4`$ of Ref. ). One observes that the distribution corresponding to fixed $`\mathrm{}`$ is very narrow whereas the Euclidean version is broad and exhibits a long tail. The striking contrast between the two histograms in Fig. 2 can be understood qualitatively. A fixed chemical distance $`\mathrm{}`$ between the origin and the trap means that there exists at least a minimum path connecting those sites whose length is $`\mathrm{}`$. Thus, the minimum time taken by a random walker to arrive at the trap site is $`t=\mathrm{}`$ independently of the lattice structure. On the other hand, a fixed Euclidean distance $`r`$ could correspond to many different chemical distances from one cluster to another. This is a consequence of the stochastic lacunarity of the fractal aggregate. A large hole between the origin and the trap implies a large minimum time to travel around the border of that hole and, obviously, a small mortality in comparison with another cluster where no holes hinder the diffusion.

Figure 3 shows the chemical mean-square displacement

$$L^2\mathrm{}^2=\underset{\text{sites}}{}\mathrm{}^2P(\mathrm{},t)$$
(26)
as a function of time. The propagator in the chemical space $`P(\mathrm{},t)`$ is obtained by summing $`P(𝐢,t)`$ over all cluster sites $`𝐢`$ on the chemical shell situated at distance $`\mathrm{}`$ from the origin. The result is compatible with the Einstein relation Eq. (11) with $`2D_{\mathrm{}}=1.20\pm 0.1`$ and $`d_w^{\mathrm{}}=2.40\pm 0.05`$. This value for $`d_w^{\mathrm{}}`$ coincides with that obtained in Ref.  and is in agreement with the value reported in Refs. .

In Fig. 4 we plot $`\mathrm{ln}(\mathrm{ln}h(\mathrm{},t))`$ versus $`\widehat{\xi }\mathrm{}/t^{1/d_w^{\mathrm{}}}`$ with $`\mathrm{}=80`$ and, according to the previous discussion, $`d_w^{\mathrm{}}=2.40`$. If the conjecture in Eq. (12) is right, we can take $`h(\mathrm{},t)\mathrm{exp}(\widehat{c}\widehat{\xi }^v)`$ as a first approximation, and hence should observe the linear behavior $`\mathrm{ln}(\mathrm{ln}h(\mathrm{},t))\mathrm{ln}\widehat{c}+v\widehat{\xi }`$ with $`\widehat{\xi }=\sqrt{2D_{\mathrm{}}}\xi `$ and $`\widehat{c}=c/(2D_{\mathrm{}})^{v/2}`$. Certainly the plot seems linear except for a portion in the range $`\widehat{\xi }2.2`$. This is a finite size effect (already analyzed in the case of the two-dimensional Sierpinski gasket in Ref. ) associated with the existence of a minimum arrival time corresponding to a random walker who travels “ballistically” along a chemical path from the origin to the trap, which in turn implies a maximum available value of $`\widehat{\xi }`$ in the simulations (in our simulations this maximum value is $`80/80^{1/d_w^{\mathrm{}}}12.9`$). A reliable interval for numerical fits should exclude this very short time regime. A linear fit in the interval $`1.6\mathrm{ln}\widehat{\xi }2.17`$, corresponding to $`200t1000`$, gives the values $`\widehat{c}=1.2\pm 0.1`$, i.e., $`c=1.3\pm 0.1`$, and $`v=1.6\pm 0.05`$. The dashed line in Fig. 4 corresponds to these values. The good agreement with numerical values in the above interval seems to assure the correctness of the approximation $`h(\mathrm{},t)h_\text{a}(\mathrm{},t)=\mathrm{exp}(\widehat{c}\widehat{\xi }^v)`$ with the values of $`c`$ and $`v`$ given above. However, the solid line in Fig. 4 is a challenge to this interpretation: one sees that the function $`h_\text{b}(\mathrm{},t)=\widehat{A}\widehat{\xi }^{\mu v}\mathrm{exp}(\widehat{c}\widehat{\xi }^v)`$ with $`v=1.70`$, $`\widehat{c}=0.9`$ (i.e. $`c=1.05`$), $`\mu =0.8`$ and $`\widehat{A}=1.1`$ is as good as $`h_\text{a}(\mathrm{},t)`$. Indeed, $`h_\text{b}(\mathrm{},t)`$ is more consistent from a theoretical point of view than $`h_\text{a}(\mathrm{},t)`$ because the expected theoretical value of $`v`$ corresponding to $`d_w^{\mathrm{}}=2.40`$ is $`v=d_w^{\mathrm{}}/(d_w^{\mathrm{}}1)=1.71`$, which is in better agreement with the exponent $`v=1.7`$ of $`h_\text{a}(\mathrm{},t)`$ than with the exponent $`v=1.6`$ of $`h_\text{b}(\mathrm{},t)`$. Finally, it should be noticed that the values $`c=1.05`$, $`v=1.70`$ are also in agreement with the corresponding parameter values of the propagator , thus supporting the guess made in Sec. II \[see below Eq. (12)\] that the dominant exponential term of the propagator and of the mortality function are the same. This leads us to consider that the set of parameters $`\widehat{c}=0.9,v=1.7,\mu =0.8`$ is more reliable than $`\widehat{c}=1.2,v=1.6,\mu =0`$. Obviously, further intensive (and extremely time consuming) computer simulations for the mortality function would be required in order to reliably determine the values of the parameters that appear in Eq. (12) and in the asymptotic corrections of $`S_N(t)`$.

We have also evaluated numerically the fractal volume in terms of the chemical distance $`V(\mathrm{})`$, i.e., the number of lattice sites inside a circumference (in chemical space) of radius $`\mathrm{}`$. The results are shown in Fig. 5. A good fit to $`dV(\mathrm{})=d_{\mathrm{}}V_0^{\mathrm{}}\mathrm{}^{d_{\mathrm{}}1}d\mathrm{}`$ is found with $`V_0^{\mathrm{}}=1.1\pm 0.2`$ and $`d_{\mathrm{}}=1.65\pm 0.05`$. Taking into account that $`d_f=91/48`$, we deduce that $`d_{\text{min}}=d_f/d_{\mathrm{}}=1.15\pm 0.05`$, which agrees with previous estimates .

## IV SIMULATION RESULTS: TERRITORY COVERED BY $`N`$ RANDOM WALKERS ON THE PERCOLATION AGGREGATE

As discussed in Sec. I, the average implicit in the evaluation of $`S_N(t)`$ is double: first, we take an average over experiments performed on the same aggregate and then a second average over different percolation clusters. The $`N`$ random walkers are always placed initially upon the site in the center of the square box. In our simulations we have performed an average over $`100`$ runs per cluster and a second average over $`2000`$ percolation clusters in order to achieve good statistics. The maximum time considered was $`t=1000`$.

According to Eqs. (13) and (14), the quotient $`S_N(t)/(\mathrm{ln}N)^\gamma `$ with $`\gamma =d_{\mathrm{}}/v`$ is only a function of $`t`$. In Fig. 6 the logarithm of that quotient is plotted versus $`\mathrm{ln}t`$ for several values of $`N`$. The data collapse and the slope close to $`0.66`$ seems to support Eq. (8) with $`\gamma =d_{\mathrm{}}/v=0.97`$, which is in agreement with similar recent results for the three-dimensional percolation aggregate . The collapse is, however, slightly poorer when the exponent $`\gamma =d_f/u=1.24`$ ($`d_f=91/48`$ and $`d_w=2.87`$ ) proposed by Havlin et al.  \[see Eq. (4)\] is used, as Fig. 6 shows. So, one migth be led to the conclusion that the correct value of $`\gamma `$ as defined above is $`d_{\mathrm{}}/v`$. But in this analysis there was no consideration of the relatively large logarithmic corrections predicted by the asymptotic analysis presented in Section II, so that the reliability of the above conclusion is seriously affected by this omission.

To illustrate this point, let us now carry out the same kind of analysis with the simulation results of $`S_N(t)`$ when the substrate is a three-dimensional Euclidean lattice. For this case it is well known  that $`S_N(t)`$ is given by an asymptotic expression with the form of Eq. (13) in which the logarithmic corrective terms are very important even for very large values of $`N`$. Indeed, the main asymptotic term leads to very poor predictions for $`S_N(t)`$, whereas the second-order approximation ($`n=2`$) gives excellent agreement with numerical simulation results. The exponent $`\gamma `$ of the main logarithmic term in $`N`$ and the time exponent $`d_{\mathrm{}}/d_w^{\mathrm{}}=d_s/2`$ are equal to $`3/2`$. We have plotted in Fig. 7 the quotient $`S_N(t)/\mathrm{ln}^\gamma N`$ versus $`\mathrm{ln}t`$ for several values of $`N`$ taking into account that the rigorous value of $`\gamma `$ is $`3/2`$. We see that the collapse is far from being perfect because the logarithmic corrections have been ignored. Nevertheless, an effective (but incorrect!) value of $`\gamma =2.75`$ yields a much better data collapse and a slope close to the theoretical value $`d_s/2=1.5`$. We thus conclude that analysis of data collapse plots based on the form of the main term of quantities such as $`S_N(t)`$ (which typically exhibit large corrective terms) should be performed with caution. The values of the exponents estimated in this way are untrustworthy because the existence of logarithmic corrections to the main term cannot simply be ignored. The value of $`\gamma =2.75`$ obtained before is then only an effective way of including all these corrective terms together but the true expression involves a main term of the form $`(t\mathrm{ln}N)^{3/2}`$ times a series similar to that given in Eq. (13). These considerations should prevent us from drawing hasty conclusions from a simple view of plots such as Figs. 6 and 7.

Finally, in Fig. 8 we show the dependence of $`S_N(t)`$ on $`N`$ and compare simulation results with the zeroth- and first-order asymptotic prediction given by Eq. (13). When the parameter set $`\widehat{c}=0.9,v=1.7,\mu =0.8,A=1`$ (see Sec. III) is used, we get results with a very familiar aspect as they are quite similar (although, perhaps the first-order approximation is too good) to that already found for Euclidean  and Sierpinski lattices . This is indeed encouraging. However, when the parameters $`\widehat{c}=1.3`$ and $`v=1.6`$ are used, we obtain a surprising and strikingly accurate zeroth-order approximation. At this point, we again suspect that this last set of parameters are only effective parameters that include the influence of the true logarithmic corrective terms in the range of $`N`$ simulated. Hence, Fig. 8 illustrates again, but from a different perspective, how the omission of important corrective terms could lead to finding effective parameters that, although providing excellent approximations in the (relatively short) range under consideration, are really erroneous.

## V Summary

In this paper, the average fractal territory covered up to time $`t`$ by $`N`$ independent random walkers all starting from the same origin on stochastic fractal lattices is calculated in terms of an asymptotic series expansion, $`_{n=0}^{\mathrm{}}_{m=0}^ns_{nm}(\mathrm{ln}N)^{d_{\mathrm{}}/vn}(\mathrm{ln}\mathrm{ln}N)^m`$ \[see Eq. (13)\], which is formally identical to those obtained for Euclidean and deterministic fractal lattices. Equation (13) is obtained by assuming that (i) the average fractal volume inside a “hypersphere” of chemical radius $`r`$ grows as $`V_0^{\mathrm{}}r^d_{\mathrm{}}`$, (ii) the distribution of the the short-time survival probability of a single random walker in the presence of a trap is narrow, so that Eq. (6) holds, and (iii) this short-time survival probability is asymptotically given by Eq. (12). We performed numerical simulations for the two-dimensional percolation aggregate at criticality which support the validity of the above assumptions. The zeroth- and first-order theoretical asymptotic expression of $`S_N(t)`$ were calculated explicitly and, for the two-dimensional percolation aggregate, they compared reasonably well with numerical simulation results. The agreement is similar to that found for non-disordered media.

In our procedure, the use of the chemical distance turns out to be fundamental because the distribution of the short-time survival probability in the chemical space is so narrow that we can safely replace the power $`N`$ of the mean value of the survival probability by the mean value of the power $`N`$ of the survival probability. (On the contrary, this does not hold at all when Euclidean distances are used.) This allowed us to easily translate the theoretical results previously derived for Euclidean and deterministic fractals  to disordered media.

###### Acknowledgements.

This work has been supported by the DGICYT (Spain) through Grant No. PB97-1501 and by the Junta de Extremadura-Fondo Social Europeo through Grant No. IPR99C031.


format: text

text:
# Self force on particle in orbit around a black hole

\[

## Abstract

We study the self force acting on a scalar charge in uniform circular motion around a Schwarzschild black hole. The analysis is based on a direct calculation of the self force via mode decomposition, and on a regularization procedure based on Ori’s mode-sum regularization prescription. We find the four self-force at arbitrary radii and angular velocities (both geodesic and non-geodesic), in particular near the black hole, where general-relativistic effects are strongest, and for fast motion. We find the radial component of the self force to be repulsive or attractive, depending on the orbit.
PACS number(s): 04.25-g, 04.70.-s, 04.70.Bw

\]

The problem of finding the equations of motion for a particle in curved spacetime has become recently extremely important, as the first generation of interferometric gravitational wave detectors will soon be operational, and with the prospects of having a gravitational wave space antenna in the not-very-distant future. The generation of very accurate templates for the waveforms detected from a system of a compact object in orbit around a supermassive black hole is an extremely hard task. It is likely that one would need to have accurate templates for as many as $`5\times 10^5`$ orbits. For such a system, accurate templates are necessary for detection, because the predicted signal-to-noise ratio for LISA is approximately of order $`10`$ for a $`1`$ year integration time. Lack of accurate templates would result in a loss of a factor of roughly the square root of the number of orbits in sensitivity , which would result in signal-to-noise ratio below the detectibility threshold.

In order to generate accurate templates, an important necessary ingredient is the inclusion of radiation reaction in the orbital evolution of the compact object. The radiation-reaction forces need to be calculated locally, i.e., in the neighborhood of the compact object. In the conventional approach, one calculates, at infinity and at the event horizon of the black hole, the fluxes of quantities, which are constants of motion in the absence of radiation reaction. Then, one uses a balance argument to relate these fluxes to the rate of change of a corresponding local quantity of the compact object. This approach generally fails because of the inadditivity of the Carter constant. For very simple cases, e.g. for circular or equatorial orbits around a Kerr black hole, the evolution of the Carter constant is trivial, such that the conventional approach is useful . However, generic orbits around a rotating black hole are neither circular nor equatorial, and consequently a new approach, which is not based on balance arguments, is of great need.

Several approaches have been suggested for the calculation of self forces. Quinn and Wald  and Mino, Sasaki, and Tanaka  recently proposed general approaches for the calculation of self forces. However, it is not presently clear how to apply these approaches directly for actual computations, the greatest problem being the calculation of the so-called “tail” term of the Green’s function, which arises from the failure of the Huygens principle in curved spacetime. (In addition, there are in general also local, Ricci-curvature coupled and Abraham-Lorentz-Dirac (ALD) type terms , which are much easier to calculate, the former vanishing identically in vacuum.)

Recently, Ori proposed a local, causal approach for the calculation of the self forces , which is based on decomposition of the self force into modes, and on a mode-sum regularization prescription (MSRP). Although MSRP is not fully proven as yet, it has already been shown to be valid for simple cases, such as scalar charges in general orbits in Schwarzschild spacetime, and, in particular, for circular orbits which we consider here. MSRP is likely to be susceptible of generalization also for massive particles in orbit around a Kerr black hole. If robust, MSRP can be of great importance for the generation of accurate templates. We hope that MSRP can be combined with other approaches, which were recently proposed, such as mode decomposition of the self forces which are sourced by just the distant past world line  or a normal-neighborhood expansion .

We first describe very briefly the main ideas of MSRP . Then, we apply MSRP for the case of a scalar particle in circular orbit around a Schwarzschild black hole, and calculate the self four-force acting on the particle linearized in its own self field.

The contribution to the physical self force from the tail part of the Green’s function can be decomposed into stationary Teukolsky modes, and then summed over the frequencies $`\omega `$ and the azimuthal numbers $`m`$. The self force equals then the limit $`ϵ0^{}`$ of the sum over all $`\mathrm{}`$ modes, of the difference between the force sourced by the entire world line (the bare force $`{}_{}{}^{\mathrm{bare}}F_{\mu }^{\mathrm{}}`$) and the force sourced by the half-infinite world line to the future of $`ϵ`$, where the particle has proper time $`\tau =0`$, and $`ϵ`$ is an event along the past ($`\tau <0`$) world line. Next, we seek a regularization function $`h_\mu ^{\mathrm{}}`$ which is independent of $`ϵ`$, such that the series $`_{\mathrm{}}({}_{}{}^{\mathrm{bare}}F_{\mu }^{\mathrm{}}h_\mu ^{\mathrm{}})`$ converges. Once such a function is found, the regularized self force is then given by $`{}_{}{}^{\mathrm{ren}}F_{\mu }^{}=_{\mathrm{}}({}_{}{}^{\mathrm{bare}}F_{\mu }^{\mathrm{}}h_\mu ^{\mathrm{}})d_\mu `$, where $`d_\mu `$ is a finite valued function. MSRP  then shows, from a local integration of the Green’s function, that the regularization function $`h_\mu ^{\mathrm{}}=a_\mu \mathrm{}+b_\mu +c_\mu \mathrm{}^1`$, and for the case of a scalar charge in circular orbit around a Schwarzschild black hole MSRP yields the values of the functions $`a_\mu ,b_\mu ,c_\mu `$ and $`d_\mu `$ analytically. In particular, it can be shown that for such orbits $`a_\mu =0=c_\mu `$ and $`d_r=0`$, such that in practice the regularization prescription of the radial force is reduced to subtracting $`b_r`$ from each $`\mathrm{}`$ mode of the bare radial force. Note that $`b_\mu `$ is just the limit $`\mathrm{}\mathrm{}`$ of $`{}_{}{}^{\mathrm{bare}}F_{\mu }^{\mathrm{}}`$.

In the following we describe the results obtained from this new approach for the case of a point-like scalar charge in circular orbit around a Schwarzschild black hole. Our results are fully relativistic, i.e., we do not introduce any simplifying assumptions such as far field or slow motion. Because of the fully relativistic nature of this study, our analysis is numerical. However, it is reasonable to expect that analytical solutions will not be available in general, except, possibly, only for very simple cases, such as static configurations . We calculate the contribution to the force which the scalar charge feels, due to its own field, to leading order in the particle’s charge. We use spherical Regge-Wheeler coordinates, for which the Schwarzschild metric is $`ds^2=\left(1\frac{2M}{r}\right)\left(dt^2+dr_{}^{}{}_{}{}^{2}\right)+r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right)`$, where $`M`$ is the black hole’s mass, and the radial Schwarzschild coordinate $`r(r^{})`$ is given implicitly by $`r^{}=r+2M\mathrm{ln}\left|r/(2M)1\right|`$. The field satisfies the wave equation $`_\mu ^\mu \mathrm{\Phi }(x^\alpha )=4\pi \rho (x^\alpha )`$, where the charge density $`\rho (x^\alpha )=q_{\mathrm{}}^{\mathrm{}}𝑑\tau \delta ^4[x^\alpha x_0^\alpha (\tau )](g)^{1/2}`$. Here, $`q`$ is the charge of the particle, whose world line is $`x_0^\alpha (\tau )`$, $`\tau `$ being its proper time, $`g`$ being the metric determinant, and $`_\mu `$ denotes covariant differentiation. We take the charge to be in circular orbit at $`r^{}=r_0^{}`$, $`\theta =\pi /2`$, and $`d\varphi /dt=\mathrm{\Omega }`$. (We are not restricted to the Keplerian angular velocity $`\mathrm{\Omega }_K`$.) We next decompose the scalar field $`\mathrm{\Phi }`$ into modes according to $`\mathrm{\Phi }=_{\mathrm{}}^{\mathrm{}}𝑑\omega _\mathrm{}me^{i\omega t}\mathrm{\Psi }_\mathrm{}m(r^{})Y^\mathrm{}m(\theta ,\varphi )/r(r^{})`$, such that the equation which governs the field $`\mathrm{\Psi }_\mathrm{}m`$ becomes

$`{\displaystyle \frac{d^2\mathrm{\Psi }_\mathrm{}m}{dr_{}^{}{}_{}{}^{2}}}`$ $`+`$ $`\left\{\omega ^2V_{\mathrm{}}[r(r^{})]\right\}\mathrm{\Psi }_\mathrm{}m=4\pi {\displaystyle \frac{q}{\gamma }}`$ (1)
$`\times `$ $`{\displaystyle \frac{\delta (rr_0)}{r_0}}\delta (\omega m\mathrm{\Omega })Y^\mathrm{}m(\pi /2,\varphi )e^{im\varphi }.`$ (2)
This equation should be solved for each mode $`\mathrm{}m`$ with boundary conditions of ingoing waves at the event horizon $`(r^{}\mathrm{})`$, and outgoing waves at infinity $`(r^{}\mathrm{})`$. The effective potential is given by $`V_{\mathrm{}}(r)=(12M/r)[2M/r^3+\mathrm{}(\mathrm{}+1)/r^2]`$, and $`\gamma =1/\sqrt{12M/rr^2\mathrm{\Omega }^2}`$. The contribution of the $`\mathrm{}m`$ mode to the force is given by $`{}_{}{}^{\mathrm{bare}}F_{\mu }^{\mathrm{}m}=q(\mathrm{\Phi }_{,\mu }^\mathrm{}m+u_\mu u^\nu \mathrm{\Phi }_{,\nu }^\mathrm{}m)`$ (note that for circular orbits $`u^\nu \mathrm{\Phi }_{,\nu }^\mathrm{}m=0`$).

We next find numerically the solutions $`\mathrm{\Psi }_\mathrm{}m^+(\mathrm{\Psi }_\mathrm{}m^{})`$ for the homogeneous equations corresponding to Eq. (2), which satisfy the boundary condition at infinity (the horizon). The components of the force are then given by

$`{}_{}{}^{\mathrm{bare}}F_{r^{}}^{\mathrm{}m}`$ $`=`$ $`2\pi q^2{\displaystyle \frac{\left|Y^\mathrm{}m(\frac{\pi }{2},0)\right|^2}{\gamma r_0^2}}\{\mathrm{Re}[W_\mathrm{}m^1(r_0^{})]\mathrm{Re}[S_\mathrm{}m(r_0^{})]`$ (3)
$`+`$ $`{\displaystyle \frac{2}{r_0}}\left(1{\displaystyle \frac{2M}{r_0}}\right)\mathrm{Re}[W_\mathrm{}m^1(r_0^{})]\mathrm{Re}[T_\mathrm{}m(r_0^{})]`$ (4)
$``$ $`{\displaystyle \frac{2}{r_0}}\left(1{\displaystyle \frac{2M}{r_0}}\right)\mathrm{Im}[W_\mathrm{}m^1(r_0^{})]\mathrm{Im}[T_\mathrm{}m(r_0^{})]`$ (5)
$`+`$ $`\mathrm{Im}[W_\mathrm{}m^1(r_0^{})]\mathrm{Im}[S_\mathrm{}m(r_0^{})]\},`$ (6)
where $`W_\mathrm{}m`$ is the Wronskian determinant of $`\mathrm{\Psi }_\mathrm{}m^{}(r^{})`$ and $`\mathrm{\Psi }_\mathrm{}m^+(r^{})`$, $`T_\mathrm{}m(r^{})=\mathrm{\Psi }_\mathrm{}m^+(r^{})\mathrm{\Psi }_\mathrm{}m^{}(r^{})`$, and $`S_\mathrm{}m(r^{})=\mathrm{\Psi }_\mathrm{}m^+(r^{})\mathrm{\Psi }_{\mathrm{}m,r^{}}^{}(r^{})+\mathrm{\Psi }_\mathrm{}m^{}(r^{})\mathrm{\Psi }_{\mathrm{}m,r^{}}^+(r^{})`$. We find that

$`{}_{}{}^{\mathrm{bare}}F_{t}^{\mathrm{}m}`$ $`=`$ $`4\pi q^2m\mathrm{\Omega }{\displaystyle \frac{\left|Y^\mathrm{}m(\frac{\pi }{2},0)\right|^2}{\gamma r_0^2}}\{\mathrm{Im}[W_\mathrm{}m^1(r_0^{})]`$ (7)
$`\times `$ $`\mathrm{Re}[T_\mathrm{}m(r_0^{})]+\mathrm{Re}[W_\mathrm{}m^1(r_0^{})]\mathrm{Im}[T_\mathrm{}m(r_0^{})]\}.`$ (8)
We also obtain $`{}_{}{}^{\mathrm{bare}}F_{\varphi }^{\mathrm{}m}=\mathrm{\Omega }^1{}_{}{}^{\mathrm{bare}}F_{t}^{\mathrm{}m}`$ and $`{}_{}{}^{\mathrm{bare}}F_{\theta }^{\mathrm{}m}=0`$. It is convenient to define new radial functions $`Z^{\pm \mathrm{}m}(r^{})`$ by $`\mathrm{\Psi }^{\pm \mathrm{}m}(r^{})=e^{\pm i\omega r^{}}Z^{\pm \mathrm{}m}(r^{})`$, which satisfy the homogeneous equations

$$\frac{d^2Z^{\pm \mathrm{}m}}{dr_{}^{}{}_{}{}^{2}}2i\omega \frac{dZ^{\pm \mathrm{}m}}{dr^{}}V_{\mathrm{}}[r^{}(r)]Z^{\pm \mathrm{}m}=0,$$
(9)
with boundary conditions $`Z^{+\mathrm{}m}(r^{}M)=1+a_1^+f_++a_2^+f_+^2+O(f_+^3)`$ and $`Z^\mathrm{}m(r^{}M)=1+a_1^{}f_{}+a_2^{}f_{}^2+O(f_{}^3)`$, where $`f_+=(\omega r)^1`$, $`f_{}=12M/r`$, and $`a_1^+=i\mathrm{}(\mathrm{}+1)/2`$, $`a_2^+=[\mathrm{}(\mathrm{}1)(\mathrm{}+1)(\mathrm{}+2)+4i\omega M]/8`$, $`a_1^{}=[1+\mathrm{}(\mathrm{}+1)]/(14i\omega M)`$, and $`a_2^{}=[\mathrm{}(\mathrm{}+1)(\mathrm{}^2+\mathrm{}+6)+4i\omega M+4]/[4(12i\omega M)(14i\omega M)]`$. We solve Eqs. (9) numerically using both Burlisch-Stoer and fourth-order Runge-Kutta integrations with adaptive step-size controls. Both integrators yield results compatible within the error limits. We place the exterior and interior boundaries at a distance of several mode wavelengths (in $`r^{}`$) from $`r_0^{}`$, and then use successive Richardson extrapolations, with increasing distance to the boundaries, until the extrapolation of the boundaries to $`r^{}\pm \mathrm{}`$ yields an error smaller than a given threshold. Notice that for modes with $`m=0`$ the wavelength is infinite, such that the boundaries cannot be taken far enough from the charge. Indeed, we find that for this case the Richardson extrapolations do not converge. Instead, we can solve for this case analytically. We find that $`f_t^{\mathrm{},m=0}=0`$ and $`f_r^{\mathrm{},m=0}=(2\pi /\gamma )(q/M)^2Y_{}^{\mathrm{}0}{}_{}{}^{2}(\pi /2,0)[2Q_{\mathrm{}}(\rho _0)dP_{\mathrm{}}(\rho _0)/d\rho +1/(1\rho _0^2)]`$, where $`\rho =(rM)/M`$. Here, $`P_{\mathrm{}},Q_{\mathrm{}}`$ are the Legendre functions of the first and second kinds, respectively. Figure 1 shows the functions $`Z_{\mathrm{}=1m=1}^\pm (r^{})`$. Similar qualitative behavior is found also for the other modes. Until the peak of the effective potential barrier the functions $`Z^{\pm \mathrm{}m}`$ vary only slowly, and then start oscillating rapidly.

The temporal component of the bare force is finite. (In the regularization scheme this corresponds to $`b_t=0`$.) MSRP predicts that $`d_t`$ exactly balances the ALD force, such that the full self force is given only by $`F_t=_{\mathrm{}}{}_{}{}^{\mathrm{bare}}F_{t}^{\mathrm{}}`$. We compare our results with their Minkowski spacetime counterparts. In flat spacetime one can solve analytically for each mode, sum over all modes, and find that $`F_t^{\mathrm{Min}}=\frac{1}{3}q^2\mathrm{\Omega }^2r_0^2\gamma _{\mathrm{Min}}^5`$. Here, $`\gamma _{\mathrm{Min}}`$ is the usual flat spacetime Lorentz factor. For the comparison we choose the same values of $`r_0,\mathrm{\Omega }`$ for the curved and flat spacetimes. Figure 2 displays $`F_t`$ as a function of $`r/M`$ for two cases: (2A) Non-geodesic circular orbits, with a fixed angular velocity $`\mathrm{\Omega }`$. (In this case the tangential velocity increases linearly with $`r`$.) When $`r/M`$ is large the value of $`F_t`$ approaches its flat spacetime counterpart. Recall that $`\mathrm{\Omega }=d\varphi /dt`$, where $`t`$ is the time of an observer at infinity. Because of the red-shift effect at small values of $`r`$, orbits with the same value of $`\mathrm{\Omega }`$ have very large proper tangential velocities at small $`r/M`$: a fixed $`\mathrm{\Omega }`$ corresponds to the ultra-relativistic limit when the orbit is close to the black hole. The second case is (2B) geodesic motion, which satisfies Kepler’s law $`\mathrm{\Omega }_K^2r^3=M`$. The innermost (unstable) causal orbit is located at $`r=3M`$. Approaching $`r=3M`$, the motion of the particle approaches the ultra-relativistic limit, which is manifested by the rapid growth of $`F_t`$. At larger radii the value of $`F_t`$ approaches the flat spacetime counterpart. This is detailed in Fig. 2(C).

Next, we study the radial, conservative component of the self force. First, we check the agreement of our numerical results with MSRP. MSRP predicts  that $`a_r=0=c_r`$ and that

$`b_r=`$ $``$ $`{\displaystyle \frac{q^2}{2r^2}}{\displaystyle \frac{1}{\gamma \sqrt{g_{tt}}}}[2{}_{2}{}^{}F_{1}^{}({\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}};1;{\displaystyle \frac{r^2\mathrm{\Omega }^2}{12M/r}})`$ (10)
$``$ $`{\displaystyle \frac{13M/r}{12M/r}}{}_{2}{}^{}F_{1}^{}({\displaystyle \frac{1}{2}},{\displaystyle \frac{3}{2}};1;{\displaystyle \frac{r^2\mathrm{\Omega }^2}{12M/r}})].`$ (11)
We check the accuracy of our numerical determination of the value of $`b_r`$ by comparison to Eq. (11). This check serves the two-fold purpose of (i) checking the numerical code, and (ii) checking the compatibility of the analytical prediction of MSRP for the regularization function $`h_r^{\mathrm{}}`$ with the numerical determination of the bare force. Figure 3(A) shows the behavior of $`{}_{}{}^{\mathrm{bare}}F_{r}^{\mathrm{}}`$, after summation over $`m`$ and $`\omega `$, as a function of $`\mathrm{}`$, and Fig. 3(B) displays $`{}_{}{}^{\mathrm{bare}}F_{r}^{\mathrm{}}b_r`$ as a function of $`\mathrm{}`$. This difference decreases like $`\mathrm{}^2`$ for large values of $`\mathrm{}`$, which confirms the predictions of MSRP for the values of $`a_r`$, $`b_r`$, and $`c_r`$. (We emphasize that we cannot test numerically the prediction of MSRP that $`d_r=0`$.) Note that the radial ALD force vanishes, such that $`{}_{}{}^{\mathrm{ren}}F_{r}^{}`$ is the full self force. Similar behavior is found also for Keplerian orbits at other values of $`r_0^{}`$, and also for non-geodesic circular motion, with angular velocities both greater or smaller than $`\mathrm{\Omega }_K`$.

The regularized component of the radial self-force $`{}_{}{}^{\mathrm{tail}}F_{r}^{\mathrm{}}`$ is obtained by subtracting $`b_r{}_{}{}^{\mathrm{bare}}F_{r}^{\mathrm{}\mathrm{}}`$ from each mode $`{}_{}{}^{\mathrm{bare}}F_{r}^{\mathrm{}}`$. The total regularized force is then obtained by

$$F_r^{\mathrm{ren}}=\underset{n=0}{\overset{\mathrm{}}{}}{}_{}{}^{\mathrm{tail}}F_{r}^{n}+_r^{\mathrm{}+1},$$
(12)
where the remainder $`_r^{\mathrm{}+1}`$ is given approximately by $`_r^{\mathrm{}+1}\mathrm{}^2({}_{}{}^{\mathrm{bare}}F_{r}^{\mathrm{}}b_r)\psi ^{(1)}(\mathrm{}+1)`$ for sufficiently large values of $`\mathrm{}`$, $`\psi ^{(1)}(x)d^2\mathrm{ln}\mathrm{\Gamma }(x)/dx^2`$ being the trigamma function. Note that for large arguments $`\psi ^{(1)}(x)x^1`$. Figure 4 displays the regularized radial self-force for Keplerian (4A) and non-Keplerian (4B) orbits. The radial self force in the far field limit ($`rM`$) is repulsive, and satisfies $`F_r^{\mathrm{ren}}\alpha q^2M^3r^5`$, where $`\alpha `$ is a dimensionless parameter of order unity. A minimum-$`\chi ^2`$ fit shows the exponent of $`r`$ to equal $`5`$ and found $`\alpha `$ to equal unity, both with $`3\%`$ errors. However, in the strong field, the force law deviates from this simple relation, and grows faster. For the non-Keplerian orbits, we find that in the slow motion limit ($`\mathrm{\Omega }\mathrm{\Omega }_K`$) the radial force is proportional to $`\mathrm{\Omega }^2`$. The exponent of $`\mathrm{\Omega }`$ is found to be $`2`$ with a $`3\%`$ error. Combined with the result for Keplerian orbits, we find that for any circular orbit, in the far field and slow motion limits, the radial force is repulsive, and is given by

$$F_r^{\mathrm{ren}}\alpha q^2(G^3/c^6)M^2\mathrm{\Omega }^2/r^2.$$
(13)
This results explains the vanishing self force in the static limit . For faster motion the $`\mathrm{\Omega }^2`$ law does not hold any more. In fact, for $`\mathrm{\Omega }>\mathrm{\Omega }_K`$ we find that the radial self force varies rapidly, and eventually changes from repulsive to attractive. The radial self force does not cause a net change in the energy of the particle. However, if the orbit has a non-zero eccentricity, this force induces an additional precession of the periastron, which in the slow motion or far field limits is retrograde. This precession has an effect on the frequencies of the emitted radiation.

Our results show that the self force can be calculated for a simple, although non-trivial, problem. MSRP was found to be useful also for other cases, e.g., static scalar and electric charges in Schwarzschild , for scalar and electric charges in circular motion in flat spacetime , and for general radial motion of scalar charges in spherical symmetry . A closely related approach was used also for a mass point in radial free fall in Schwarzschild . We hope that similar methods can be used for more realistic cases, which may be relevant for the orbital evolution of compact objects around black holes.

I thank Scott Hughes and Kip Thorne for valuable discussions. I am indebted to Amos Ori for many stimulating discussions and for letting me use his results before their publication. This research was supported by NSF grants AST-9731698 and PHY-9900776, and by NASA grant NAG5-6840.


text:
# 1 Introduction

## 1 Introduction

In spite of the fact that the top quark has been discovered already several years ago its interactions are still unknown. Therefore it remains an open question if the top-quark couplings obey the Standard Model (SM) scheme of the electroweak forces or there exists a contribution from physics beyond the SM. In this talk♯1♯1♯1Presented by B. Grzadkowski I will try to use angular and energy distributions of top-quark decay products in the process $`e^+e^{}t\overline{t}\mathrm{}^\pm /\stackrel{()}{b}\mathrm{}`$ in order to estimate how precisely top quark couplings could be determined at future linear collider.

We will parameterize $`t\overline{t}`$ couplings to the photon and the $`Z`$ boson in the following way

$$\mathrm{\Gamma }_{vt\overline{t}}^\mu =\frac{g}{2}\overline{u}(p_t)\left[\gamma ^\mu \{A_v+\delta A_v(B_v+\delta B_v)\gamma _5\}+\frac{(p_tp_{\overline{t}})^\mu }{2m_t}(\delta C_v\delta D_v\gamma _5)\right]v(p_{\overline{t}}),$$
(1)
where $`g`$ denotes the $`SU(2)`$ gauge coupling constant, $`v=\gamma ,Z`$, and

$$A_\gamma =\frac{4}{3}\mathrm{sin}\theta _W,B_\gamma =0,A_Z=\frac{1}{2\mathrm{cos}\theta _W}\left(1\frac{8}{3}\mathrm{sin}^2\theta _W\right),B_Z=\frac{1}{2\mathrm{cos}\theta _W}$$

denote the SM contributions to the vertices. Among the above non-SM form factors, $`\delta A_{\gamma ,Z}`$, $`\delta B_{\gamma ,Z}`$, $`\delta C_{\gamma ,Z}`$ describe $`CP`$-conserving while $`\delta D_{\gamma ,Z}`$ parameterizes $`CP`$-violating interactions.

Similarly, we will adopt the following parameterization of the $`Wtb`$ vertex suitable for the $`t`$ and $`\overline{t}`$ decays:

$`\mathrm{\Gamma }_{Wtb}^\mu ={\displaystyle \frac{g}{\sqrt{2}}}V_{tb}\overline{u}(p_b)\left[\gamma ^\mu (f_1^LP_L+f_1^RP_R){\displaystyle \frac{i\sigma ^{\mu \nu }k_\nu }{M_W}}(f_2^LP_L+f_2^RP_R)\right]u(p_t),`$
$`\overline{\mathrm{\Gamma }}_{Wtb}^\mu ={\displaystyle \frac{g}{\sqrt{2}}}V_{tb}^{}\overline{v}(p_{\overline{t}})\left[\gamma ^\mu (\overline{f}_1^LP_L+\overline{f}_1^RP_R){\displaystyle \frac{i\sigma ^{\mu \nu }k_\nu }{M_W}}(\overline{f}_2^LP_L+\overline{f}_2^RP_R)\right]v(p_{\overline{b}}),`$
(2)
where $`P_{L/R}=(1\gamma _5)/2`$, $`V_{tb}`$ is the $`(tb)`$ element of the Kobayashi-Maskawa matrix and $`k`$ is the momentum of $`W`$.

It will be assumed here that interactions of leptons with gauge bosons are properly described by the SM. Through the calculations all fermions except the top quark will be considered as massless. We will also neglect terms quadratic in non-standard form factors.

## 2 Angular and Energy Distributions

In this section we will present results for $`d^2\sigma /dx_fd\mathrm{cos}\theta _f`$ of the top-quark decay product $`f`$, where $`f`$ could be either $`\mathrm{}^\pm `$ or $`\stackrel{()}{b}`$, $`x_f`$ denotes the normalized energy of $`f`$ and $`\theta _f`$ is the angle between the $`e^{}`$ beam direction and the direction of $`f`$ momentum in the $`e^+e^{}`$ CM frame.

Using the technique developed by Kawasaki, Shirafuji and Tsai and adopting the general formula for the $`t\overline{t}`$ distribution $`d\sigma (s_+,s_{})/d\mathrm{\Omega }_t`$ found by Brzezinski et al., one obtains the following result for the distribution:

$$\frac{d^2\sigma }{dx_fd\mathrm{cos}\theta _f}=\frac{3\pi \beta \alpha _{\text{EM}}^2}{2s}B_f\left[\mathrm{\Theta }_0^f(x_f)+\mathrm{cos}\theta _f\mathrm{\Theta }_1^f(x_f)+\mathrm{cos}^2\theta _f\mathrm{\Theta }_2^f(x_f)\right],$$
(3)
where $`\beta `$ is the top velocity, $`\alpha _{\text{EM}}`$ is the fine structure constant and $`B_f`$ denotes the appropriate branching fraction. The energy dependence is specified by the functions $`\mathrm{\Theta }_i^f(x_f)`$, explicit forms of which could be found in a recent paper by the authors of this contribution. They are parameterized both by production and decay form factors.

The angular distribution $`d\sigma /d\mathrm{cos}\theta _f`$ for $`f`$ could be easy obtained from eq.(3) by the integration over the energy of $`f`$. It turns out that for $`f=\mathrm{}`$ all the non-standard contributions from the decay vertex disappear upon integration over the energy $`x_f`$.

The fact that the angular leptonic distribution is insensitive to corrections to the $`V`$-$`A`$ structure of the decay vertex allows for much more clear tests of the production vertices through a measurement of the distribution, since that way we can avoid a contamination from non-standard structure of the decay vertex. As an illustration, we define a $`CP`$-violating asymmetry which could be constructed using the angular distributions of $`f`$ and $`\overline{f}`$:

$$𝒜_{CP}(\theta _f)=\left[\frac{d\sigma ^+(\theta _f)}{d\mathrm{cos}\theta _f}\frac{d\sigma ^{}(\pi \theta _f)}{d\mathrm{cos}\theta _f}\right]/\left[\frac{d\sigma ^+(\theta _f)}{d\mathrm{cos}\theta _f}+\frac{d\sigma ^{}(\pi \theta _f)}{d\mathrm{cos}\theta _f}\right],$$
(4)
where $`d\sigma ^{+/}`$ is referring to $`f`$ and $`\overline{f}`$ distributions, respectively. The asymmetry for $`f=\mathrm{}`$ is sensitive to $`CP`$ violation originating exclusively from the production mechanism, i.e. $`CP`$-violating form factors $`\delta D_\gamma `$ and $`\delta D_Z`$ while the decay vertex enters with the SM $`CP`$-conserving coupling. For bottom quarks the effect of the modification of the decay vertex is contained in corrections to so called depolarization factors, $`\alpha ^b+\alpha ^{\overline{b}}\mathrm{Re}(f_2^R\overline{f}_2^L)`$ with SM $`CP`$-conserving contribution from the production process. As it is seen from fig.1 for $`\sqrt{s}=1\text{TeV}`$ the asymmetry could be quite large, e.g., reaching for the semileptonic decays $`20\%`$ for $`\mathrm{Re}(\delta D_\gamma )=\mathrm{Re}(\delta D_Z)=0.2`$.

## 3 Optimal Observable Analysis

Using the double angular and energy distributions an expected statistical uncertainty for determination of real parts for all the form factors has been found adopting optimal observables and varying the beam polarizations $`P_e^{}`$ and $`P_{e^+}`$. $`|\mathrm{cos}\theta |0.9`$ has been assumed as a cut for the polar angle. For $`t\overline{t}`$ tagging efficiency in $`\mathrm{}`$ \+ 4 jet channel we adopted 60% and for the integrated luminosity we chose the TESLA design with $`L=500\text{fb}^1`$ at $`\sqrt{s}=500\text{GeV}`$.

Generically we have observed that positive polarization led to smaller statistical errors for the eight form factors in the production vertices. For each form factor we have adjusted the optimal beam polarization such that the statistical error was minimal:

$$\begin{array}{cc}\mathrm{\Delta }[\text{Re}(\delta A_\gamma )]=0.16\hfill & \mathrm{for}P_e^{}=0.7\mathrm{and}P_{e^+}=0.7\hfill \\ \mathrm{\Delta }[\text{Re}(\delta A_Z)]=0.07\hfill & \mathrm{for}P_e^{}=0.5\mathrm{and}P_{e^+}=0.4\hfill \\ \mathrm{\Delta }[\text{Re}(\delta B_\gamma )]=0.09\hfill & \mathrm{for}P_e^{}=0.2\mathrm{and}P_{e^+}=0.2\hfill \\ \mathrm{\Delta }[\text{Re}(\delta B_Z)]=0.27\hfill & \mathrm{for}P_e^{}=0.4\mathrm{and}P_{e^+}=0.4\hfill \\ \mathrm{\Delta }[\text{Re}(\delta C_\gamma )]=0.11\hfill & \mathrm{for}P_e^{}=0.1\mathrm{and}P_{e^+}=0.0\hfill \\ \mathrm{\Delta }[\text{Re}(\delta C_Z)]=1.11\hfill & \mathrm{for}P_e^{}=0.1\mathrm{and}P_{e^+}=0.0\hfill \\ \mathrm{\Delta }[\text{Re}(\delta D_\gamma )]=0.08\hfill & \mathrm{for}P_e^{}=0.2\mathrm{and}P_{e^+}=0.1\hfill \\ \mathrm{\Delta }[\text{Re}(\delta D_Z)]=14.4\hfill & \mathrm{for}P_e^{}=0.2\mathrm{and}P_{e^+}=0.1\hfill \end{array}$$
(5)
As it is seen the precision of $`\delta \{C,D\}_Z`$ measurement would be very poor even for the optimal polarization. In addition, determination of $`\delta D_\gamma `$ would be difficult as well since its error varies rapidly with the polarization. For example, $`\mathrm{\Delta }[\text{Re}(\delta D_\gamma )]`$ becomes 0.86 for $`P_e^{}=0.1/P_{e^+}=0.1`$ and 0.99 for $`P_e^{}=0.3/P_{e^+}=0.1`$. The source of that sensitivity is hidden in the neutral current structure with $`\mathrm{sin}^2\theta _W0.23`$. Indeed, the optimal polarization becomes $`P_e^{}=0.1`$ instead of 0.2 $`(\mathrm{\Delta }[\text{Re}(\delta D_\gamma )]=0.09)`$ for $`\mathrm{sin}^2\theta _W=0.25`$. On the other hand, a good determination (almost independently of the polarization) could be expected for $`f_2^R`$. Indeed, the best precision is

$$\mathrm{\Delta }[\text{Re}(f_2^R)]=0.01\mathrm{for}P_e^{}=0.8\mathrm{and}P_{e^+}=0.8$$
(6)
however, we have $`\mathrm{\Delta }[\text{Re}(f_2^R)]=0.03`$ even for $`P_e^{}=P_{e^+}=0`$.

## 4 Summary and Conclusions

We have presented here the angular and energy distributions for $`\stackrel{()}{f}`$ in the process $`e^+e^{}t\overline{t}\stackrel{()}{f}\mathrm{}`$, where $`f=\mathrm{}`$ or $`b`$ quark. The most general ($`CP`$-violating and $`CP`$-conserving) couplings for $`\gamma t\overline{t}`$, $`Zt\overline{t}`$ and $`Wtb`$ have been assumed. The bottom-quark mass has been neglected and we have kept only terms linear in anomalous couplings.

Test of CP violation has also been discussed, introducing a CP-sensitive asymmetry $`𝒜_{CP}`$ as an example.

Using the double angular and energy distribution of a lepton we have found that at $`\sqrt{s}=500\text{GeV}`$ with the integrated luminosity $`L=500\text{fb}^1`$ the best determined top-quark coupling would be the axial coupling of the Z boson with the error $`\mathrm{\Delta }[\text{Re}(\delta A_Z)]=0.07`$ while the lowest precision is expected for $`\text{Re}(\delta D_Z)`$ with $`\mathrm{\Delta }[\text{Re}(\delta D_Z)]=14.4`$.

## Note added

After this work has been presented, a paper by Rindani appeared where the double angular and energy leptonic distribution has been also found.

## Acknowledgments

BG is grateful to the organizers of the PASCOS99 conference for creating a very warm and inspiring atmosphere during the meeting. This work is supported in part by the State Committee for Scientific Research (Poland) under grant 2 P03B 014 14 and by Maria Skłodowska-Curie Joint Fund II (Poland-USA) under grant MEN/NSF-96-252.


text:
# Untitled Document

A simple question about a complicated object

The complicated object is the cohomological induction functor (for which the biblical reference is Knapp & Vogan \[KV\]). Here is the simple question. For definiteness sake set $`G:=PU(n,1)`$, $`H:=PU(k,1)`$ with $`kn`$, and let $`_\rho ^n`$ (resp. $`_\rho ^k`$) be the category of Harish-Chandra modules with the generalized infinitesimal character of the trivial module attached to $`G`$ (resp. $`H`$). I’ll surprise nobody by claiming that $`H`$ is a subgroup of $`G`$. What’s far less obvious, but proved by Khoroshkin \[K\], is the existence of a full embedding $`F`$ of $`_\rho ^k`$ into $`_\rho ^n`$ (Fuser \[F,Thm I.4.2\] showed that $`F`$ is even Ext-full, that is compatible with Ext calculus), prompting the question : is there a geometric interpretation of the embedding $`F`$ ? The first candidate for $`F`$ is the (ordinary) induction functor ; but this fails miserably — so let’s break the pseudo-suspense of this introduction by saying that I claim that $`F`$ is (isomorphic to) a certain cohomological induction functor, and conjecture that this phenomenon is general.

1. Statements

Let $`G`$ be a center free connected semisimple Lie group, $`KG`$ a maximal compact subgroup, $`𝔤𝔨`$ the respective complexified Lie algebras. Let’s start by recalling the notion of Harish-Chandra module. Say that a $`𝔤`$-module $`V`$ is $`𝔨`$-finite if it is a sum of finite dimensional sub-$`𝔨`$-modules, and that $`V`$ is an $`(𝔤,𝔨)`$-module if it is $`𝔨`$-finite and $`𝔨`$-semisimple. The category $`=(𝔤,K)`$ of Harish-Chandra modules is the full subcategory of $`𝔤`$-mod whose objects are those $`(𝔤,𝔨)`$-modules of finite length $`V`$ such that for any finite dimensional $`𝔨`$-invariant subspace $`FV`$ the action of $`𝔨`$ on $`F`$ exponentiates to $`K`$. The category $``$ is a $`\mathrm{}`$-category in the sense of Bass \[B\] page 57. Let $`I`$ be the annihilator of the trivial module in the center of $`U(𝔤)`$, let

$$_\rho =_\rho (𝔤,K)$$

be the full sub-$`\mathrm{}`$-category of $``$ whose objects are annihilated by some power of $`I`$, let $``$ be the set of isomorphism classes of simple objects of $`_\rho `$ \[it is a finite set\] ; for each $`i`$ choose a representative $`V_ii`$ and let $`\mathrm{}(i)`$ be the projective dimension of $`V_i`$ \[i.e. the supremum in $`\mathrm{}\{+\mathrm{}\}`$ of the set $`\{n\mathrm{}|Ext^n(V_i,)0\}`$\].

(1) Definition. The $`_\rho `$-ordering is the smallest partial ordering $``$ on $``$ satisfying

$$\begin{array}{c}i,j\\ \\ \mathrm{}(j)=\mathrm{}(i)+1<\mathrm{}\\ \\ Ext^1(V_j,V_i)0\end{array}\}ij.$$

(2) Definition. The sub-$`\mathrm{}`$-category generated by the subset $`𝒥`$ of $`_\rho `$ is the full sub-$`\mathrm{}`$-category $`𝒥__\rho `$ of $`_\rho `$ characterized by the condition that an object $`V`$ of $`_\rho `$ belongs to $`𝒥__\rho `$ iff each simple subquotient of $`V`$ is isomorphic to $`V_j`$ for some $`j𝒥`$.

(3) Definition. Say that a full sub-$`\mathrm{}`$-category $`𝒞`$ of $`_\rho `$ is Ext-full in $`_\rho `$ if for all $`V,W𝒞`$ the natural morphism

$$Ext_𝒞^{}(V,W)Ext__\rho ^{}(V,W)$$

is an isomorphism.

For $`i`$ put $`𝒥_i:=\{j|ji\}`$, let $`\theta `$ be the Cartan involution of $`(𝔤,K)`$, denote by $`d`$ the dimension of $`G/K`$, and consider the following

(4) Property of $`G`$. For each $`i`$ such that $`V_i`$ is unitary the cohomology $`H^{d\mathrm{}(i)}(𝔤,K;V_i)`$ is nonzero and there is a $`\theta `$-stable parabolic subalgebra of $`𝔤`$ with Levi subgroup $`L=L_i`$ (see Vogan \[V2,4.1,4.2\] for definitions) satisfying

(a) the corresponding cohomological induction functor $`F`$ (see \[KV\]) sets up an equivalence

$$_\rho (𝔩,LK)\stackrel{}{}𝒥_i_{_\rho (𝔤,K)};$$

(b) $`F\mathrm{}V_i;`$

(c) $`𝒥_i_{_\rho (𝔤,K)}`$ is Ext-full in $`_\rho (𝔤,K)`$ ;

(d) if $`a`$ is nonzero vector of $`H^{d\mathrm{}(i)}(𝔤,K;V_i)`$ and $`V`$ a simple object of $`_\rho (𝔩,LK)`$, then the map

$$\begin{array}{ccc}H^{}(𝔩,LK;V)& & H^{d\mathrm{}(i)+}(𝔤,K;FV)\\ \\ x& & F(x)a\end{array}$$

\[where $``$ denotes the cup-product\] is an isomorphism \[of $`H^{}(𝔩,LK;\mathrm{})`$-modules\] ;

(e) we have $`2\mathrm{}(i)=d+dimL/(LK)`$.

Note once and for all that (e) follows from (a) by the well known argument which consists in setting $`V:=\mathrm{}`$ and using Poincaré duality.

(6) Conjecture. All center free connected semisimple Lie groups have Property (4).

A partial proof (with explicitly indicated gaps) of the fact that $`PU(n,1)`$, $`P\text{Spin}(n,1)`$ and $`SL(3,\mathrm{})`$ have Property (4) is contained in the expanded version of this text, downloadable from

http://www.iecn.u-nancy.fr/$``$gaillard/Recherche/Ci/ci.html

* * *

\[F\] Fuser A., Autour de la conjecture d’Alexandru, Thèse de l’Univesité Nancy 1 (1997).

\[K\] Khoroshkin S.M., Category of Harish-Chandra modules of the group SU(n,1). Funct. Anal. Appl. 14 (1980) 153-155.

\[KV\] Knapp A. & Vogan D., Cohomological induction and unitary representations, Princeton University Press (1995).

\[V\] Vogan D., Cohomology and group representations, Proc. Symp. Pure Math. 61 (1997) 219-243.


text:
# The Energy Spectrum of TeV Gamma-Rays from the Crab Nebula as measured by the HEGRA system of imaging air Čerenkov telescopes

## 1 Introduction

The Crab Nebula has been observed and studied over an enormously broad photon energy range embracing the radio, optical, and X-ray bands, as well as, high energy $`\gamma `$-ray region up to hundreds of TeV. The various theoretical scenarios of photon emission are primarily based on the Synchro Compton model (Gould 1965) which combines the synchrotron and inverse Compton (IC) emissions from high-energy electrons, which are accelerated up to $`100`$ TeV and interact with the magnetic field and the low frequency seed photons within the Nebula (de Jager & Harding 1992; Atoyan & Aharonian 1996; de Jager et al. 1996; Hillas et al. 1998). The predicted IC spectrum in the TeV energy domain appears to be very sensitive to the model parameters, such as the value of the magnetic field, the nature of the seed photons, the maximum energy of electrons, etc. Although IC scenarios of the photon emission are widely believed to be appropriate for the Crab Nebula one can not exclude the possible contribution of $`\gamma `$-ray fluxes from $`\pi ^{}`$-decay (see Atoyan & Aharonian 1996, Bednarek & Protheroe 1997). Thus, a measurement of TeV Crab Nebula spectrum sets major constraints on theoretical expectations and precise spectral measurements allow to fix the model parameters.

The imaging air Čerenkov technique was successfully used for observations of the Crab Nebula in TeV $`\gamma `$-rays. Since the time of detection at a 9$`\sigma `$ confidence level by the Whipple group (Weekes et al. 1989) a number of observations of the Crab Nebula have been made at TeV energies (for reviews see Ong 1998; Catanese & Weekes 1999). By now the Crab Nebula is established as the standard candle of steady TeV $`\gamma `$-ray emission. The optical Nebula of the Crab has an angular extension of about 6 arc min (Hester et al. 1995). The standard analysis for the HEGRA system of IACTs provides an angular resolution of about $`0.1^{}`$ (see Figure 1). In the present analysis the Crab Nebula was assumed to be as a $`\gamma `$-ray point source (see Figure 2). The detailed mapping of the TeV $`\gamma `$-rays from the Crab Nebula undertaken by Hofmann (1999) provided an upper limit for the Crab angular extension in TeV $`\gamma `$-rays of $`1.5`$ arcmin. No pulsed emission has been seen so far from the Crab (Gillanders et al. 1997; Burdett et al. 1999; Aharonian et al. 1999a). Vacanti et al (1991) provided a first measurement of the energy spectrum of TeV $`\gamma `$-ray emission from the Crab Nebula using 10 m Whipple telescope data and derived a power-law spectrum index of 2.4 and a differential $`\gamma `$-ray flux at 400 GeV of $`2.510^{10}\mathrm{photons}\mathrm{cm}^2\mathrm{s}^1`$. Note that the absolute calibration of individual imaging air Čerenkov telescopes is quite difficult because there is no “test” beam of TeV $`\gamma `$-rays. The uncertainties in absolute calibration propagate to uncertainties of estimated $`\gamma `$-ray fluxes. Different methods of telescope calibration have been developed recently in order to reduce systematic errors which could influence the estimate of the absolute $`\gamma `$-ray flux and the slope of energy spectrum (see e.g., Frass et al. 1997). The energy spectrum measurements heavily rely on Monte Carlo simulations of the telescope response. Thus in the past the analysis of the same observational data using various methods of telescope calibration revealed, not infrequently, very different estimates of the telescope energy threshold and of the TeV $`\gamma `$-ray fluxes. These estimates may vary by more than a factor of 2. Thus, based on data taken with the prototype HEGRA imaging Čerenkov telescope (CT2), Konopelko et al. (1996) detected the TeV $`\gamma `$-ray signal from the Crab Nebula with a signal-to-noise ratio of 10$`\sigma `$ and derived from the data a power-law index of the energy spectrum of $`2.7\pm 0.1`$ and an integral flux of TeV $`\gamma `$-rays above 1 TeV of $`8(\pm 1)_{\mathrm{Stat}}(\pm 2.4)_{\mathrm{Syst}}10^{12}\mathrm{photons}\mathrm{cm}^2\mathrm{s}^1`$. However, later on, the detailed treatment of the telescope hardware (optical smearing, photon-to-photoelectrons conversion efficiency, etc) allowed a more precise estimate of the TeV $`\gamma `$-ray flux using the single telescope data, as $`1.5(\pm 0.2)_{\mathrm{Stat}}(+1.00.5)_{\mathrm{Syst}}10^{11}\mathrm{photons}\mathrm{cm}^2\mathrm{s}^1`$ above 1 TeV (Petry et al. 1996).

Significant improvements in the telescope hardware as well as in the simulations and the data analysis (Fegan, 1997) recently provided measurements of the Crab Nebula energy spectrum over the energy range from 200 GeV up to 50 TeV by several groups, Whipple (Hillas et al. 1998), CANGAROO (Tanimori et al. 1998), and CAT (see Barrau 1998) exploiting the imaging Čerenkov technique. An additional measurement was provided by a low energy air shower array (Amenomori et al. 1999). However the uncertainties of the $`\gamma `$-ray flux estimates as well as of the energy spectrum slope remain rather large, in particular at the high energy end of the measurements i.e. beyond 10 TeV. For example the Tibet data (Amenomori et al. 1999) show significantly higher (by a factor of 2) $`\gamma `$-ray fluxes compared with the results obtained using the IACTs in the energy range from 3 to $`18`$ TeV, and are in favor of a gradual steepening of the spectral slope at high energies. Thus additional precise measurements of the Crab Nebula energy spectrum in TeV $`\gamma `$-rays are of great importance. The HEGRA system of IACT provides such data. It was primarily designed for detailed spectral measurements in the TeV energy domain utilizing the advantages of stereoscopic observations. Stereo imaging gives several advantages for spectral studies, compared to a single telescope: (i) direct measurement of the shower impact parameter with an accuracy better than 10 m (ii) good energy resolution of 18% (iii) wide dynamic range from 500 GeV to 20 TeV (iv) extended abilities for systematic studies using several images for an individual shower. The detailed systematic studies for the spectrum evaluation technique have been recently made using Mrk 501 1997 observational data (Aharonian et al. 1999b, 1999c). The performance of the system was discussed by Konopelko et al. (1999a). Here we present Crab Nebula data taken with the HEGRA telescope system in 1997/1998 and 1998/1999 observational campaigns. The data were analyzed using a new technique of energy spectrum evaluation for the stereoscopic observations. We also discuss the physics implications of the present results for the modeling of the TeV $`\gamma `$-ray emission from the Crab Nebula.

## 2 Observational data

The Crab Nebula was extensively observed with the HEGRA IACT system in two observing seasons from September 1997 to March 1998 and from October 1998 to April 1999. The observations were made with the stereoscopic system of IACTs which are located on La Palma, Canary Islands (Aharonian et al. 1999c). Each of the telescopes consists of a 8.5 $`\mathrm{m}^2`$ reflector focussing the Čerenkov light onto a photomultiplier tube camera. The 271 photomultipliers in the camera were arranged in a hexagonal matrix covering a field of view with a radius of $`2.15^{}`$. The telescope camera was triggered when the signal in two next neighbors of the 271 photo multiplier tubes exceeded a threshold of 8 photoelectrons, and the system readout started when at least two telescopes were triggered by Čerenkov light from an air shower. The detection rate was 12.6 Hz near the Zenith in December 1997 and dropped down to about 10 Hz in December 1999 for the 4-telescope system due to aging of the PMTs and reduced mirror reflectivity.

The Crab Nebula was observed in a “wobble mode”; i.e., the telescopes were pointed in Declination $`\pm 0.5^{}`$ away from the nominal Crab Nebula position (the sign of the angular shift was altered from one run of 20 min to the next). This is useful for continuous monitoring of the cosmic-ray background because it positions the OFF-source region symmetric to the camera center, and $`1^{}`$ apart from the ON-source region. Observations of the Crab Nebula at zenith angles up to 50 degree were made with 4 telescopes from 1997 September 1 to 1998 March 29, for a total of 82.5 hr of data taken at good weather. Through the fire at the HEGRA site one of the telescopes was damaged and was out of operation for a month in October-November 1997. At that time, Crab Nebula observations were made with only three telescopes in the system, providing an event rate of 10.3 Hz near the zenith. Due to unstable weather and a substantial amount of dust which came from Sahara desert to the island, the average detection rate in 1998 February-March was reduced down to 10.7 Hz in observations near the zenith. At the beginning of 1998/1999 observational period 4 telescopes were operational. Since October 1998 the HEGRA collaboration operates 5 telescopes. However, for technical reasons one telescope (CT2) was out of operation since December 1998 until the end of the first observation campaign (April 1999). During the last observational period (since August 1999) observations at zenith angles less than $`50^{}`$ were taken for about 76.1 hrs. In addition, the observations at large zenith angles (LZA) ($`50^{}<\theta <65^{}`$), were carried out for a total of 24 hrs in order to study the performance of the telescope system at LZA and to extend the measurements of the Crab Nebula energy spectrum beyond 10 TeV. The total exposure times for the three periods are summarized in Table 1.

Only data taken under good weather conditions were used in the analysis. In order to exclude data taken under less than optimal telescope performance conditions the entire database has been checked very carefully as follows: First, each night the compressed protocols of the system performance were transfered to one of the collaboration host institutes where they were scanned by software tools which closely monitor the status of the telescopes’ hardware (single pixel rates, trigger rates of system telescopes, tracking accuracy, etc). This information was accumulated in a corresponding database which was used afterwards for a standard data reduction procedure. The final condensed data file for each particular run contains all the information needed for data analysis. In addition a specific software tool was developed which allows to control a posteriori for each data run (i) the system trigger rate, taking into account the zenith angle dependence (ii) the angular shape of the cosmic ray images, tested by a $`\chi ^2`$–criterion for the deviation of the mean scaled Width distribution for a single run from the corresponding average distribution filled over an extended sample of runs (iii) the flatness of the $`\theta ^2`$–distribution for the isotropic cosmic ray images over the full field of view (iv) the image Size distributions for each individual telescope.

## 3 Analysis

The stereoscopic imaging analysis of the data is based on the geometrical reconstruction of the shower arrival direction and the shower core position in the observation plane, as well as on the joint parameterization of the shape of the Čerenkov light images. The simultaneous registration of several ($``$2) Čerenkov light images from an air shower provides an angular resolution of $`0.^{}1`$ for $`\gamma `$-ray showers. For each individual shower, stereoscopic observations allow to determine the position of the shower axis. Thus, at first only air showers within a certain impact distance $`\mathrm{R}_0`$ from the center of the telescope system were selected. The limiting upper radius of $`\mathrm{R}_0=200\mathrm{m}`$ was used for zenith angles less than 50 degrees, and a significantly larger radius of 400 m for the large zenith angle observations ($`>50`$ degrees). The effective collection area in observations at LZA dramatically increases at high energies, far beyond the limiting radius of 200 m (Konopelko et al. 1999b). For the data taken at zenith angles up to $`50^{}`$ an orientation cut $`\theta ^2<0.05[\mathrm{deg}^2]`$ was applied, where $`\theta ^2`$ is the squared angular distance of the reconstructed source position from the true source position. In addition the data were analyzed using the mean scaled Width parameter, $`<\stackrel{~}{w}>`$. To compensate for the dependence of the image shape on primary shower energy and distance from shower core to the telescope (impact parameter), the standard parameter Width (Fegan, 1997) ($`w^k`$), calculated for each telescope, is scaled according to the Monte Carlo predicted values, for $`\gamma `$-rays $`<w>_{ij}^k`$, taken for the corresponding bin of reconstructed distance from the telescope to the shower core (i) and for the corresponding bin of image size (total number of photoelectrons in the image)(j) (Aharonian et al. 1999b, 1999c). The mean scaled Width parameter is defined for each individual shower as follows

$$<\stackrel{~}{w}>=1/N\underset{k=1}{\overset{N}{}}w^k/<w>_{ij}^k$$
(1)
where $`N`$ is the number of triggered telescopes. This parameter was introduced in order to provide an almost constant $`\gamma `$-ray acceptance over the dynamic energy range of the telescope system. The optimum cut on mean scaled Width is about 1.1, which gives a $`\gamma `$-ray acceptance of $`60`$% at Small Zenith Angle (SZA). However, for a precise determination of $`\gamma `$-ray spectra, a loose cut on mean scaled Width ($`<\stackrel{~}{w}><1.2`$) has been so far used in the data analysis in order to maximize the $`\gamma `$-ray acceptance and to minimize systematic errors related to cut efficiencies. Thus, the second $`\gamma `$-ray selection criterion was $`<\stackrel{~}{w}>1.2`$. This set of cuts was found to be optimal for spectrum studies (Aharonian et al. 1999b, 1999c). These loose analysis cuts provide a Crab Nebula $`\gamma `$-ray rate of 83 $`\gamma `$s/hr at SZA (less that $`25^{}`$) for the 5-IACT system. The corresponding energy threshold of the $`\gamma `$-rays is about 500 GeV. For the LZA data the looser orientation cut of $`\theta ^2<0.1[\mathrm{deg}^2]`$ was used because of the lower accuracy of the arrival direction reconstruction for the $`\gamma `$-ray showers. In order to improve the cosmic ray rejection in observations at LZA an additional parameter, mean scaled Length, $`<\stackrel{~}{l}>`$ was used, which is defined by analogy with $`<\stackrel{~}{w}>`$ (Konopelko et al. 1999b). These two parameters, $`<\stackrel{~}{w}>`$ and $`<\stackrel{~}{l}>`$, can be used for calculating a Mahalanobis distance, MD (Mahalanobis 1963), in two-dimensional space as

$$\mathrm{MD}=((1<\stackrel{~}{\mathrm{w}}>)^2/\sigma _{<\stackrel{~}{\mathrm{w}}>}^2+(1<\stackrel{~}{\mathrm{l}}>)^2/\sigma _{<\stackrel{~}{\mathrm{l}}>}^2)^{1/2}$$
(2)
where $`\sigma _{<\stackrel{~}{w}>}`$ and $`\sigma _{<\stackrel{~}{l}>}`$ are the standard deviations for the corresponding distributions of $`<\stackrel{~}{w}>`$ and $`<\stackrel{~}{l}>`$. The optimum value of the MD cut for LZA is found to be 1.5. Note that this analysis improves the enhancement factor by $`30`$% (it gives $`50`$% acceptance of $`\gamma `$-rays) in observations at LZA, whereas it gives only marginal improvement for the data taken at SZA. The Crab Nebula $`\gamma `$-ray rate in observations at LZA ($`60^{}`$) is about 16 $`\gamma `$’s/hr with a corresponding energy threshold of $``$5 TeV. Note that SZA observations give a $`\gamma `$-ray rate at high energies (above 3 TeV) of $`8\gamma `$s/hr. A summary of the data is shown in Table 2.

The observations of the Crab Nebula have been made during 6 periods which differ in the system configuration, reflectivity of the mirrors, light reflection by the pixel funnels and camera protecting plate etc. All that affects the hardware event rate of the telescope system, $`\mathrm{R}_{\mathrm{exp}}`$. These changes of system performance were implemented in the Monte Carlo simulations. Assuming the standard chemical composition of the primary cosmic rays (Wiebel, 1994) the calculated detection rates, $`\mathrm{R}_{\mathrm{MC}}`$, were adjusted to the measured rates (see Table 3).

The collection areas, as a function of energy and zenith angle, for $`\gamma `$-ray showers have been inferred from Monte Carlo simulations (Konopelko et al. 1999a). The rms error of the energy determination is $`\mathrm{\Delta }\mathrm{E}/\mathrm{E}0.18`$. The Monte Carlo studies show that for a good energy resolution of 18% this approach does not distort the initial spectrum shape. The collection area for $`\gamma `$-rays rises very quickly in the energy range near the energy threshold of the telescope system, which is 500 GeV, whereas it is almost constant at the energies above $`3`$ TeV. Even slight variations of the trigger threshold could lead to noticeable systematic changes in the predicted spectral behavior in the energy range of $`0.51`$ TeV. This effect leads to a noticeable probability for “sub-threshold” triggers. In addition, the trigger level for different camera pixels is slightly different even after very accurate adjustment of the high voltage using the calibration laser runs. Measurements of the trigger setting for a number of camera pixels revealed variation in the trigger threshold of order $`10`$%. These variations were implemented into the simulations in order to estimate the corresponding systematic error of the energy spectrum at energies below 1 TeV. The fine tuning of the Monte Carlo simulations with respect to the IACT system data provided measurements of the flux of the cosmic ray protons in the energy range from 1.3 to 10 TeV (Aharonian et al. 1999d). The proton fluxes as measured by the HEGRA IACT system are perfectly consistent with the results of a bulk of satellite experiments held in this energy range.

The procedure for the evaluation of the energy spectrum using the stereoscopic observations was discussed in detail by Aharonian et al. (1997); Hofmann (1997), and, more recently in Aharonian et al. (1999c). In the stereoscopic observations the impact distance of the shower axis to a system telescope can be measured with an accuracy $`10`$ m. The energy E of a $`\gamma `$-ray shower is defined by interpolation over the “size” parameter $`\mathrm{S}`$ (total number of photoelectrons in Čerenkov light image) at a fixed impact distance $`R`$, as $`\mathrm{E}=\mathrm{f}_{\mathrm{MC}}(\mathrm{S},\mathrm{R},\theta )`$, where $`\theta `$ is the zenith angle and $`\mathrm{f}_{\mathrm{MC}}`$ is a function obtained from Monte Carlo simulations. Note that the Monte Carlo simulations used here include the sampling of detector response in great detail (Hemberger 1998). The energy distribution for the ON- and OFF-source events, after the orientation and shape image cuts, were histogrammed over the energy range from 500 GeV to 30 TeV with 8 bins per decade. The $`\gamma `$-ray energy spectrum was obtained by subtracting ON- and OFF-histograms and dividing the resulting energy distribution by the corresponding collection area and the $`\gamma `$-ray acceptance. In the present Crab Nebula analysis the energy spectrum measurements were extended up to large zenith angles ($`65^{}`$). The data were processed independently for each of the four zenith angle bins: $`(0^{}25^{}),(25^{}40^{}),(40^{}50^{}),(50^{}65^{})`$. The corresponding effective collection areas as well as the cut efficiencies were calculated as a function of the zenith angle. First, the energy spectra were derived for all zenith angle bins independently. Note that the spectra evaluated at different zenith angles are in a good agreement. For the final energy spectrum the different zenith angle bins were joined according to the prescription:

$`dJ_\gamma ^i/dE={\displaystyle \underset{j=1}{\overset{4}{}}}w_j(dJ_\gamma ^i/dE)_j\mathrm{\Theta }(E^iE_{th}^j),`$
$`w_j=t_j/t_0,i=1,n;`$ (3)
where $`dJ_\gamma ^i/dE`$, $`(dJ_\gamma ^i)_j/dE`$ are the differential energy spectra at energy $`E^i`$ as measured over all zenith angle ranges, and for the particular zenith angle bin (j), respectively. $`E_{th}^j`$ is an estimated energy threshold for the zenith angle bin $`j`$, $`t_j`$ is the observation time for the $`j`$-bin on the zenith angle, and $`t_0`$ is the total observation time. The first three zenith angle bins were joined using the time dependent weights $`\mathrm{w}_\mathrm{j}=\mathrm{t}_\mathrm{j}/\mathrm{t}_0`$. Finally the spectrum measured in the zenith angle range of $`(0^{}50^{})`$ was combined with the spectrum derived from large zenith angle data, $`(50^{}65^{})`$ using the weights based on the estimated statistical errors for both spectra. Such procedure takes into account the advantageous $`\gamma `$-ray rate in LZA observations at high energies.

The statistics of the $`\gamma `$-rays from the Crab Nebula provides a measurement of the energy spectrum up to a few tens of TeV. However, detection of Čerenkov light images with extremely large amplitudes - several thousands of ph.e. - is complicated by the nonlinearity in the PMT response as well as by the saturation in the 8 bit Flash-ADC readout. Measurements of the photomultiplier response under high light loads over the extended sample of the EMI 9073 PMTs gave a calibration function which was used to correct the image amplitudes. The readout of the HEGRA IACT is based on the sampling of Čerenkov light time impulse by the 16 FADC bins of $`8`$ ns each (Hess et al. 1998). The time pulses from the air showers with a full width at half maximum of a few ns were widened using an electronic scheme in order to fit into several FADC bins for the accurate measurement of the time profile. The smoothing of the FADC signal was unfolded back to the impulse, which almost always fits 2 FADC bins. The calibrated amplitude, summed over two FADC bins, is used as a measure of the pixel signal. For the high energy air showers the FADC signals run into saturation and the simple unfolding procedure fails. For such pulses the initial amplitude is reconstructed using the additional calibration function obtained by simultaneous measurements of light flashes with FADCs and a 14 bit ADC. This procedure drastically extends the dynamic range of the FADC readout.

To avoid the saturation problem one might only use images detected from air showers at large impact distances from the telescope system (e.g. beyond 150 m). The size of these images is very small even for high energy events because of the low Čerenkov light density far off the shower axis. However these images are very often truncated by the camera edge and do not allow a proper reconstruction of the shower impact point and of the shower energy. This effect becomes less important in observations at LZA because of the high shower maximum height (the images shrink to the camera center). In the present analysis the maximum impact distance of the shower core from the center of the system was extended up to 400 m for observations at LZA. Observations at LZAs permit measurements of the energy spectrum far beyond 10 TeV. The images of $`\gamma `$-ray air showers observed at LZAs have small Size and are not influenced by the saturation effect.

## 4 Results

We have observed the Crab Nebula extensively in two observational seasons with the HEGRA IACT system. The HEGRA system of 5 IACTs currently has a sensitivity which allows the detection of a $``$5$`\sigma `$ signal from the Crab Nebula within 1 hr of observation time (see Figure 1). The integral $`\gamma `$-ray fluxes measured during the different observational periods are consistent within the estimated statistical and systematic error (see Table 3). The differential energy spectrum of the Crab Nebula has been derived from the HEGRA data for two observational campaigns using recently developed advanced techniques for the measurements of the spectrum using stereoscopic data taken at small and large zenith angles. The Crab Nebula differential energy spectrum derived from SZA data matches quite well the spectrum derived at LZA (see Figure 3). The $`\gamma `$-ray rate measured at energies above 10 TeV in observations at LZAs exceeds the corresponding rate measured at SZA by a factor of 3. The LZA data are not affected by saturation effects. At 3.7$`\sigma `$ confidence level 27 $`\gamma `$-ray events from the Crab Nebula were detected in the highest energy bin from 17.8 to 23.7 TeV (see Figure 4). One may expect that a number of $`\gamma `$-ray events in the highest energy bin are spilled over from the lower energies. However, given the good energy resolution of 20% and the power law energy spectrum, such effect is very small and is compensated almost by the backwards influx of the $`\gamma `$-rays from the energies above (see e.g., Aharonian et al. 1995). Finally, the simulations show that spilling over of low energy $`\gamma `$-rays does not influence the resulting fluxes measured at the upper end of the power law spectrum.

The analysis for the different system configurations as well as for different trigger threshold values gives a differential energy spectrum of the Crab Nebula measured at zenith angles up to $`60^{}`$

$$\mathrm{dJ}_\gamma /\mathrm{dE}=(2.79\pm 0.02\pm 0.5)10^7(\frac{\mathrm{E}}{1\mathrm{TeV}})^{2.59\pm 0.03\pm 0.05}\mathrm{ph}\mathrm{m}^2\mathrm{s}^1\mathrm{TeV}^1$$
(4)
The statistical and systematic errors are also given. The final Crab Nebula spectrum as measured by the HEGRA collaborations is shown in Figure 4. The measured $`\gamma `$-ray fluxes are given in Table 4. The Crab Nebula energy spectrum is best fitted by a pure power law in the energy range 1-20 TeV. It does not exclude a possible slight steepening of the energy spectral usually predicted by inverse Compton modeling of TeV $`\gamma `$-ray emission. A fit with a logarithmic steepening of the power law spectrum gives the following result

$$\mathrm{dJ}_\gamma /\mathrm{dE}=(2.67\pm 0.01\pm 0.5)10^7(\frac{\mathrm{E}}{1\mathrm{TeV}})^{2.47\pm 0.1\pm 0.05(0.11\pm 0.10)\mathrm{log}(\mathrm{E})}\mathrm{ph}\mathrm{m}^2\mathrm{s}^1\mathrm{TeV}^1$$
(5)
Such a fit indicates the slight flattening of the spectrum at low energies as predicted by the IC calculations. However, the change of the energy spectrum slope is within the current statistical and systematic errors, and the data for the overall energy range are consistent with a simple power law fit in 0.5-20 TeV. The HEGRA Crab Nebula data match well the recent 20 TeV data published by the CANGAROO group (Tanimori et al. 1998), and are consistent with a flat power law index of $`2.5`$ beyond 20 TeV. The compilation of the world data is given in Figure 5. All data are consistent within statistical and systematic errors, except possibly for the Tibet data which show relatively higher fluxes.

## 5 Astrophysics implications

The TeV energy spectrum of the Crab Nebula as measured by the HEGRA system of IACTs is consistent with the expectations for the TeV $`\gamma `$-ray emission from pulsar-driven Nebulae (plerions). According to this scenario the ultra relativistic electrons, accelerated in the pulsar wind shock, produce TeV $`\gamma `$-rays through the IC scattering with soft photons within the Nebula. The predicted fluxes of TeV $`\gamma `$-rays rely on the spatial distribution of the magnetic field within the Nebula as determined by the parameter $`\sigma `$ (ratio of energy density of magnetic field to the particle energy density) and/or by the average magnetic field $`<\mathrm{B}>`$ in the optical nebula. The HEGRA data are shown in Figure 6 together with predicted spectra using two SSC models of TeV $`\gamma `$-ray emission (de Jager et al. 1996; Atoyan & Aharonian 1996). One may conclude that both models fit the HEGRA data rather well.

According to the calculations of the TeV $`\gamma `$-ray emission by de Jager et al. (1996) the $`\gamma `$-ray flux from the Crab Nebula at TeV energies constrains the choice of the parameter $`\sigma `$. The IC spectrum computed by de Jager et al. (1996), assuming for the parameter $`\sigma `$ a value $`0.003`$ gives a good fit to the HEGRA data. This value of the parameter $`\sigma `$ corresponds to the best-fitting magneto-hydrodynamic (MHD) solution of electron propagation in the Crab Nebula as found by Kennel & Coroniti (1984).

In another approach, assuming the spatial distribution of magnetic field in the Crab Nebula, one can determine the average magnetic field $`<\mathrm{B}>`$ in the optical nebula (Gould 1965). According to the calculations of Atoyan & Aharonian (1996), made within the framework of the MHD model of Kennel & Coroniti (1984), the average magnetic field $`<\mathrm{B}>`$ is determined by the TeV $`\gamma `$-ray flux as $`<\mathrm{B}>\mathrm{J}_\gamma ^{0.5}10^5\mathrm{G}`$ where $`\mathrm{J}_\gamma ,\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1\mathrm{TeV}^1`$ is a differential flux at 1 TeV. Thus the HEGRA spectrum gives an average magnetic field strength $`<\mathrm{B}>(1.7\pm 0.3)10^4\mathrm{G}`$. This value is consistent with the estimate derived by Hillas et al. (1998) from the Crab Nebula date taken with the 10 m Whipple telescope.

The IC energy spectrum of $`\gamma `$ rays from the Crab Nebula, measured in the energy range from 1 TeV to 10 TeV, is likely to be a power law $`\mathrm{dJ}_\gamma /\mathrm{dE}\mathrm{E}^\alpha =\mathrm{E}^{2.6}`$. At the same time, for the energy range above 10 TeV calculations predict gradual steepening with $`\alpha `$2.7 and 2.9 at 10 and 30 TeV, respectively. That is due to both the energy loss of the ultra high energy electrons by fast synchrotron cooling and the Klein-Nishina effect in the cross-section of the inverse Compton scattering. Atoyan & Aharonian (1996) and Bednarek & Protheroe (1997) have shown that $`\pi ^0`$-decay $`\gamma `$-ray fluxes, due to the relativistic protons accelerated in the Crab Nebula, may noticeably contribute at energies above 10 TeV. However the HEGRA Crab Nebula data expanded up to 20 TeV are still consistent with the pure IC spectrum. To assess the contribution of $`\pi ^0`$-produced $`\gamma `$-rays from the Crab Nebula measurements above 30 TeV are needed. Note that the LZA technique could help to perform such observations.

The predicted IC $`\gamma `$-ray spectrum of the Crab Nebula is rather flat in the energy range below $`1`$ TeV. It could be well approximated by $`\mathrm{dJ}_\gamma /\mathrm{dE}\mathrm{E}^{2.0}`$ at 100 GeV. Detection of a gradual flattening in this energy range will prove the SSC scenario of TeV $`\gamma `$-ray emission. However the low energy points at the HEGRA Crab Nebula spectra ($`\mathrm{E}_\gamma <1\mathrm{TeV}`$) are strongly affected by possible systematic errors ($`50`$%) and do not allow such a conclusion. Future observations of the Crab Nebula with the forthcoming low threshold high sensitivity Čerenkov detectors (see Catanese & Weekes 1999) will offer precise measurements in these energy range.

## Acknowledgments

The support of the German ministry for Research and technology BMBF and of the Spanish Research Council CYCIT is gratefully acknowledged. We thank the Instituto de Astrophysica de Canarias for the use of the site and for supplying excellent working conditions at La Palma. We gratefully acknowledge the technical support staff of the Heidelberg, Kiel, Munich, and Yerevan Institutes.