text:
# The Transition from Inspiral to Plunge for a Compact Body in a Circular Equatorial Orbit Around a Massive, Spinning Black Hole

## I INTRODUCTION AND SUMMARY

The space-based Laser Interferometer Space Antenna (LISA) , if it flies, is likely to detect and study the gravitational waves from white dwarfs, neutron stars and small black holes with masses $`\mu 1M_{}`$, spiraling into Massive ($`M10^5`$$`10^8M_{}\mu `$) black holes in the nuclei of distant galaxies . In preparation for these studies, it is necessary to understand, theoretically, the radiation-reaction-induced evolution of the inspiral orbits, and the gravitational waveforms that they emit.

Regardless of an orbit’s shape and orientation, when $`\mu M`$ the orbital evolution can be divided into three regimes: (i) The adiabatic inspiral regime, in which the body gradually descends through a sequence of geodesic orbits with gradually changing “constants” of the motion $`E=(`$energy), $`L=(`$polar component of angular momentum), and $`Q=(`$Carter constant). (ii) A transition regime, in which the character of the orbit gradually changes from inspiral to plunge. (iii) A plunge regime, in which the body plunges into the horizon along a geodesic with (nearly) unchanging $`E`$, $`L`$ and $`Q`$.

The plunge regime, being (essentially) ordinary geodesic motion, is well understood; and the adiabatic inspiral regime is the focus of extensive current research (see, e.g. ). By contrast, so far as we are aware, there have been no publications dealing with the transition regime.

In 1990-91, we carried out an initial exploration of the transition regime for the special case of circular, equatorial orbits; but due to the press of other projects we did not publish it. Now, with prospects for LISA looking good, we have resurrected our work and present it here, in the context of LISA.

We begin, in Sec. II, by summarizing some key, well-known details of the inspiral and plunge regimes. Then in Sec. III A we present a qualitative picture of the transition from inspiral to plunge, based on the motion of a particle in a slowly changing effective potential (Fig. 1). With the aid of this qualitative picture, in Sec. III B we derive a non-geodesic equation of motion for the transition regime, and in Sec. III C we construct the solution to that equation of motion (Figs. 2 and 3). Then in Sec. IV, with the aid of our solution, we estimate the gravitational-wave signal strength from the transition regime and the signal-to-noise ratio that it would produce in LISA. We conclude that, with some luck, LISA may be able to detect and study the transition waves. In Sec. V we make concluding remarks about the need for further research.

## II ADIABATIC INSPIRAL AND PLUNGE

Throughout this paper we use Boyer-Lindquist coordinates $`(t,r,\theta ,\varphi )`$  for the massive hole’s Kerr metric, and we use geometrized units, with $`G=c=1`$. The hole’s mass is $`M`$ and the inspiraling body’s mass is $`\mu \eta M`$. We use $`M`$ and $`\mu `$ to construct dimensionless versions (denoted by tildes) of many dimensionfull quantities; for example, $`\stackrel{~}{r}=r/M`$, and $`\stackrel{~}{t}=t/M`$. The hole’s dimensionless spin parameter is $`a(`$spin angular momentum$`)/M^2`$ (with $`1<a<+1`$). The body moves around its circular, equatorial orbit in the $`+\varphi `$ direction, so $`a>0`$ corresponds to an orbit that is prograde relative to the hole’s spin, and $`a<0`$ to a retrograde orbit.

When the inspiraling body is not too close to the innermost stable circular orbit (isco), it moves on a circular geodesic orbit with dimensionless angular velocity

$$\stackrel{~}{\mathrm{\Omega }}M\mathrm{\Omega }=\frac{d\varphi }{d\stackrel{~}{t}}=\frac{1}{\stackrel{~}{r}^{3/2}+a}$$
(1)
(where $`\varphi `$ is angle around the orbit) and with orbital energy

$$E=\eta M\frac{12/\stackrel{~}{r}+a/\stackrel{~}{r}^{3/2}}{\sqrt{13/\stackrel{~}{r}+2a/\stackrel{~}{r}^{3/2}}}.$$
(2)
As it moves, the body radiates energy into gravitational waves at a rate given by

$$\dot{E}_{\mathrm{GW}}=\dot{E}=\frac{32}{5}\eta ^2\stackrel{~}{\mathrm{\Omega }}^{10/3}\dot{},$$
(3)
where $`\dot{}`$ is a general relativistic correction to the Newtonian, quadrupole-moment formula (Table II of Ref. ). This energy loss causes the orbit to shrink adiabatically at a rate given by

$$\frac{dr}{dt}=\frac{\dot{E}_{\mathrm{GW}}}{dE/dr}.$$
(4)
The inspiral continues adiabatically until the body nears the isco, which is at the dimensionless radius $`\stackrel{~}{r}_{\mathrm{isco}}=r_{\mathrm{isco}}/M`$ given by

$`\stackrel{~}{r}_{\mathrm{isco}}=`$ $`3+Z_2\mathrm{sign}(a)[(3Z_1)(3+Z_1+2Z_2)]^{1/2},`$ (7)
$`Z_11+(1a^2)^{1/3}[(1+a)^{1/3}+(1a)^{1/3}],`$
$`Z_2(3a^2+Z_{1}^{}{}_{}{}^{2})^{1/2};`$
cf. Table I. The circular geodesic orbit at the isco has dimensionless angular velocity (Table I), energy, and angular momentum given by

$$\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}M\mathrm{\Omega }=\frac{1}{\stackrel{~}{r}_{\mathrm{isco}}^{}{}_{}{}^{3/2}+a},$$
(8)
$$\stackrel{~}{E}_{\mathrm{isco}}\frac{E_{\mathrm{isco}}}{\mu }=\frac{E_{\mathrm{isco}}}{\eta M}=\frac{12/\stackrel{~}{r}_{\mathrm{isco}}+a/\stackrel{~}{r}_{\mathrm{isco}}^{}{}_{}{}^{3/2}}{\sqrt{13/\stackrel{~}{r}_{\mathrm{isco}}+2a/\stackrel{~}{r}_{\mathrm{isco}}^{}{}_{}{}^{3/2}}},$$
(9)
$$\stackrel{~}{L}_{\mathrm{isco}}\frac{L_{\mathrm{isco}}}{\mu M}=\frac{L_{\mathrm{isco}}}{\eta M^2}=\frac{2}{\sqrt{3\stackrel{~}{r}_{\mathrm{isco}}}}\left(3\sqrt{\stackrel{~}{r}_{\mathrm{isco}}}2a\right).$$
(10)
As the body nears the isco, its inspiral gradually ceases to be adiabatic and it enters the transition regime (Sec. III). Radiation reaction (as controlled by $`\dot{E}_{\mathrm{GW}}`$) continues to drive the orbital evolution throughout the transition regime, but gradually becomes unimportant as the transition ends and pure plunge takes over.

The plunge is described to high accuracy by reaction-free geodesic motion; Eqs. (33.32) of Ref. . Up to fractional corrections of order $`\eta ^{4/5}`$, the orbital energy and angular momentum of the plunging body are equal to $`E_{\mathrm{isco}}`$ and $`L_{\mathrm{isco}}`$ throughout the plunge. \[cf. Eq. (40) below\].

## III THE TRANSITION FROM ADIABATIC INSPIRAL TO PLUNGE

### A Qualitative Explanation of Transition

As the body nears its innermost stable circular orbit, $`r=r_{\mathrm{isco}}`$, the adiabatic approximation begins to break down. This breakdown can be understood in terms of the effective potential, which governs geodesic radial motion via the equation

$$\left(\frac{d\stackrel{~}{r}}{d\stackrel{~}{\tau }}\right)^2=\left(\frac{dr}{d\tau }\right)^2=\stackrel{~}{E}^2V(\stackrel{~}{r},\stackrel{~}{E},\stackrel{~}{L}),$$
(11)
where $`\stackrel{~}{E}E/\mu =E/(\eta M)`$, $`\stackrel{~}{L}L/(\mu M)=L/(\eta M^2)`$, and $`\stackrel{~}{\tau }\tau /M`$ are the body’s dimensionless energy, angular momentum, and proper time. The explicit form of the effective potential can be inferred from Eqs. (33.32) and (33.33) of MTW:

$`V(\stackrel{~}{r},\stackrel{~}{E},\stackrel{~}{L})`$ $`=`$ $`\stackrel{~}{E}^2{\displaystyle \frac{1}{\stackrel{~}{r}^4}}([\stackrel{~}{E}(\stackrel{~}{r}^2+a^2)\stackrel{~}{L}a]^2`$ (13)
$`(\stackrel{~}{r}^22\stackrel{~}{r}+a^2)[\stackrel{~}{r}^2+(\stackrel{~}{L}\stackrel{~}{E}a)^2]).`$
For a Schwarzschild black hole, this reduces to

$$V(\stackrel{~}{r},\stackrel{~}{E},\stackrel{~}{L})=\left(1\frac{2}{\stackrel{~}{r}}\right)\left(1+\frac{\stackrel{~}{L}^2}{\stackrel{~}{r}^2}\right)\text{for }a=0$$
(14)
(cf. Eq. (25.16) of MTW ).

Throughout the inspiral and transition regimes, the body moves along a nearly circular orbit; its change of radius during each circuit around the hole is $`\mathrm{\Delta }rr`$. (Only after the body is well into its final plunge toward the hole does $`\mathrm{\Delta }r`$ become comparable to $`r`$.) This near-circular motion guarantees that the ratio of the energy radiated to angular momentum radiated is equal to the body’s orbital angular velocity :

$$\frac{d\stackrel{~}{E}}{d\stackrel{~}{\tau }}=\stackrel{~}{\mathrm{\Omega }}\frac{d\stackrel{~}{L}}{d\stackrel{~}{\tau }}.$$
(15)
Correspondingly, in and near the transition regime, which occupies a narrow range of radii around $`\stackrel{~}{r}_{\mathrm{isco}}`$, the body’s energy and angular momentum are related by*** In reality, finite-mass-ratio effects, including those discussed in the paragraph preceding Eq. (22) below, will alter these energy-angular-momentum relations by amounts that scale as the first and higher powers of $`\eta `$. For example, in going from Eq. (15) to (16), there can be an integration constant $`\delta \stackrel{~}{E}`$ (which scales as $`\eta `$ or some higher power) so $`\stackrel{~}{E}=\stackrel{~}{E}_{\mathrm{isco}}+\delta \stackrel{~}{E}+\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}\xi `$. In the presence of such effects, we redefine $`\stackrel{~}{r}_{\mathrm{isco}}`$, $`\stackrel{~}{E}_{\mathrm{isco}}`$, and $`\stackrel{~}{L}_{\mathrm{isco}}`$ to be the values of these parameters at which the $`\eta `$-corrected $`V(\stackrel{~}{r},\stackrel{~}{E},\stackrel{~}{L})`$ has a flat inflection point, as in Fig. 1, and $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}`$ to be the orbital angular velocity at this $`\stackrel{~}{r}_{\mathrm{isco}}`$. Then Eqs. (16) remain valid even for finite mass ratio $`\eta `$.

$$\stackrel{~}{E}=\stackrel{~}{E}_{\mathrm{isco}}+\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}\xi ,\stackrel{~}{L}=\stackrel{~}{L}_{\mathrm{isco}}+\xi .$$
(16)
By combining Eqs. (16) and (13), we can regard the body’s effective potential as a function of $`\stackrel{~}{r}`$ and the difference $`\xi \stackrel{~}{L}\stackrel{~}{L}_{\mathrm{isco}}`$ of its orbital angular momentum from that of the isco.

Figure 1 shows $`V(\stackrel{~}{r},\xi )`$ for a sequence of angular momenta $`\xi _1,\mathrm{},\xi _5`$ around $`\xi =0`$. As $`\xi `$ decreases to $`\xi =0`$, the minimum of the potential flattens out and disappears; and just when it is disappearing, the minimum’s radius $`r_{\mathrm{min}}`$ is moving inward at an infinite rate: $`dr_{\mathrm{min}}/d\xi \mathrm{}`$ as $`\xi 0`$.

In the adiabatic regime of large $`\xi `$, the body sits always at the minimum of the effective potential. Its orbit is a slowly shrinking circle, guided inward by the motion of the minimum. As $`\xi `$ nears zero and the minimum’s inward speed grows large, the body’s inertia prevents it from continuing to follow the minimum. The body begins to lag behind, as depicted at $`\xi =\xi _2`$ in Fig. 1. This lag invalidates the adiabatic inspiral analysis of Sec. II and initiates the transition regime.

As $`\xi `$ continues to decrease, there comes a point (near $`\xi _5`$ in Fig. 1) at which the effective potential has become so steep that its inward force on the body dominates strongly over radiation reaction. There the transition regime ends, and the body begins to plunge inward rapidly on a nearly geodesic orbit with nearly constant $`\stackrel{~}{E}`$ and $`\stackrel{~}{L}`$. The objectives of the following subsections are to derive a set of equations describing the transition regime (Sec. III B), and show how the transition matches smoothly onto the adiabatic regime at large positive $`\xi `$ and to the plunge regime at large negative $`\xi `$ (Sec. III C).

### B Equation of Motion for Transition Regime

Throughout the transition regime, because the body moves on a nearly circular orbit with radius close to $`r_{\mathrm{isco}}`$, and because the body’s small mass $`\mu \eta MM`$ keeps its radiation reaction weak, its angular velocity remains very close to $`\mathrm{\Omega }_{\mathrm{isco}}`$

$$\frac{d\varphi }{d\stackrel{~}{t}}\stackrel{~}{\mathrm{\Omega }}\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}},$$
(17)
and its proper time ticks at very nearly the standard isco rate

$$\frac{d\stackrel{~}{\tau }}{d\stackrel{~}{t}}\left(\frac{d\stackrel{~}{\tau }}{d\stackrel{~}{t}}\right)_{\mathrm{isco}}=\frac{\sqrt{13/\stackrel{~}{r}_{\mathrm{isco}}+2a/\stackrel{~}{r}_{\mathrm{isco}}^{}{}_{}{}^{3/2}}}{1+a/\stackrel{~}{r}_{\mathrm{isco}}^{}{}_{}{}^{3/2}};$$
(18)
cf. Eq. (5.4.5a) of .

This nearly circular motion at $`\stackrel{~}{r}\stackrel{~}{r}_{\mathrm{isco}}`$ produces gravitational waves which carry off angular momentum and energy at very nearly the same rate as they would for circular geodesic motion at $`\stackrel{~}{r}_{\mathrm{isco}}`$. This means that $`\stackrel{~}{E}`$ and $`\stackrel{~}{L}`$ evolve in accord with Eqs. (16), where

$$\frac{d\xi }{d\stackrel{~}{\tau }}=\kappa \eta ,$$
(19)
and

$$\kappa =\frac{32}{5}\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}^{}{}_{}{}^{7/3}\frac{1+a/\stackrel{~}{r}_{\mathrm{isco}}^{}{}_{}{}^{3/2}}{\sqrt{13/\stackrel{~}{r}_{\mathrm{isco}}+2a/\stackrel{~}{r}_{\mathrm{isco}}^{}{}_{}{}^{3/2}}}\dot{}_{\mathrm{isco}};$$
(20)
cf. Eqs. (3), (16), (18), and Table I. It is the smallness of $`\eta \mu /M`$ (e.g., $`\eta =10^5`$ for the realistic case of a $`10M_{}`$ black hole spiraling into the $`10^6M_{}`$ black hole) that makes the angular momentum $`\xi `$ evolve very slowly and keeps the body in a nearly circular orbit throughout the transition regime \[cf. the factors of $`\eta `$ that appear in Eqs. (19) and (33)—which with Eqs. (22) and (35) imply $`d\stackrel{~}{r}/d\stackrel{~}{\tau }\eta ^{3/5}`$.\]

In the transition regime, the body’s radial motion is described by the geodesic equation of motion with a radial self force per unit massThis radial self force, like the radiation reaction force that drives the inspiral, is produced by interaction of the body with its own gravitational field—that field having been influenced by the black hole’s spacetime geometry; see, e.g., Ref. . The contravariant radial component of the self force, with dimensionality restored using $`r=M\stackrel{~}{r}`$ and $`\tau =M\stackrel{~}{\tau }`$, is $`(dp^r/d\tau )_{\mathrm{self}}=(\mu d^2r/d\tau ^2)_{\mathrm{self}}=(\mu /M)(d^2\stackrel{~}{r}/d\stackrel{~}{\tau }^2)_{\mathrm{self}}=\eta ^2\stackrel{~}{F}_{\mathrm{self}}`$.  $`\eta \stackrel{~}{F}_{\mathrm{self}}`$ inserted on the right-hand side:

$$\frac{d^2\stackrel{~}{r}}{d\stackrel{~}{\tau }^2}=\frac{1}{2}\frac{V(\stackrel{~}{r},\xi )}{\stackrel{~}{r}}+\eta \stackrel{~}{F}_{\mathrm{self}}.$$
(21)
(We write it as $`\eta \stackrel{~}{F}_{\mathrm{self}}`$ because its magnitude is proportional to $`\eta =\mu /M`$.)

The radial self force $`\eta \stackrel{~}{F}_{\mathrm{self}}`$ is nondissipative (since it has hardly any radial velocity with which to couple). This contrasts with the $`\varphi `$-directed radiation-reaction force, which couples to the orbital angular velocity to produce a shrinkage of the body’s angular momentum \[Eq. (19)\] and a corresponding decrease of its energy, $`d\stackrel{~}{E}/d\stackrel{~}{\tau }=\stackrel{~}{\mathrm{\Omega }}d\xi /d\stackrel{~}{\tau }`$. Because the radial force is nondissipative, it is of little importance. It can be absorbed into the nondissipative effective potential term $`\frac{1}{2}V/\stackrel{~}{r}`$ in the equation of motion. Doing so will not change the general character of the effective potential, as depicted in Fig. 1; it will merely change, by fractional amounts proportional to $`\eta `$, the various parameters that characterize the effective potential: the location $`\stackrel{~}{r}_{\mathrm{isco}}`$ of the innermost stable circular orbit (at which the $`\xi =0`$ effective potential curve has its inflection point), the values at the isco of the orbital energy and angular momentum $`\stackrel{~}{E}_{\mathrm{isco}}`$ and $`\stackrel{~}{L}_{\mathrm{isco}}`$, and the constant $`\alpha `$ defined below. There will be other $`\mathrm{O}(\eta )`$ changes in $`\stackrel{~}{r}_{\mathrm{isco}}`$, $`\stackrel{~}{E}_{\mathrm{isco}}`$, $`\stackrel{~}{L}_{\mathrm{isco}}`$ and $`\alpha `$ caused by the body’s own perturbation of the hole’s spacetime geometry. In this paper, we shall ignore all such changes, and correspondingly we shall neglect the radial self force $`\eta \stackrel{~}{F}_{\mathrm{self}}`$.

We shall describe the body’s location in the transition regime by

$$R\stackrel{~}{r}\stackrel{~}{r}_{\mathrm{isco}}.$$
(22)
Throughout the transition regime both $`R`$ and $`\xi `$ are small, and correspondingly the effective potential can be expanded in powers of $`R`$ and $`\xi `$. Up through cubic terms in $`R`$ and linear terms in $`\xi `$ (the order needed for our analysis), the effective potential takes the form

$$V(R,\xi )=\frac{2\alpha }{3}R^32\beta R\xi +\text{constant},$$
(23)
where $`\alpha `$ and $`\beta `$ are positive constants that we shall evaluate below. Note that for $`\xi =0`$, this is a simple cubic potential with inflection point at $`R=0`$, i.e. at $`\stackrel{~}{r}=\stackrel{~}{r}_{\mathrm{isco}}`$; and note that for $`\xi >0`$, it acquires a maximum and a minimum, while for $`\xi <0`$ it is monotonic; cf. Fig. 1. By inserting Eq. (23) into Eq. (21), setting $`\stackrel{~}{r}=\stackrel{~}{r}_{\mathrm{isco}}+R`$, and neglecting the radial self force or absorbing it into $`\stackrel{~}{r}_{\mathrm{isco}}`$, $`\alpha `$ and $`\beta `$ as described above, we obtain the following radial equation of motion:

$$\frac{d^2R}{d\stackrel{~}{\tau }^2}=\alpha R^2+\beta \xi .$$
(24)
By then setting $`\stackrel{~}{\tau }0`$ at the moment when $`\xi =0`$ and using Eq. (19) for the rate of change of $`\xi `$, so

$$\xi =\eta \kappa \stackrel{~}{\tau },$$
(25)
we bring our equation of motion into the form

$$\frac{d^2R}{d\stackrel{~}{\tau }^2}=\alpha R^2\eta \beta \kappa \stackrel{~}{\tau }.$$
(26)
We shall explore the consequences of this equation of motion in the next subsection, but first we shall deduce the values of $`\alpha `$ and $`\beta `$.

The constants $`\alpha `$ and $`\beta `$ can be evaluated from the following relations, which follow directly from (23), (22) and (16):

$$\alpha =\frac{1}{4}\left(\frac{^3V(\stackrel{~}{r},\stackrel{~}{E},\stackrel{~}{L})}{\stackrel{~}{r}^3}\right)_{\mathrm{isco}},$$
(27)
$$\beta =\frac{1}{2}\left(\frac{^2V(\stackrel{~}{r},\stackrel{~}{E},\stackrel{~}{L})}{\stackrel{~}{L}\stackrel{~}{r}}+\stackrel{~}{\mathrm{\Omega }}\frac{^2V(\stackrel{~}{r},\stackrel{~}{E},\stackrel{~}{L})}{\stackrel{~}{E}\stackrel{~}{r}}\right)_{\mathrm{isco}}.$$
(28)
By inserting expression (13) into these relations, we obtain $`\alpha `$ and $`\beta `$ in the limit $`\eta \mu /M0`$:

$`\alpha `$ $`=`$ $`{\displaystyle \frac{3}{\stackrel{~}{r}_{\mathrm{isco}}^6}}\left(\stackrel{~}{r}^2+2[a^2(\stackrel{~}{E}^21)\stackrel{~}{L}^2]\stackrel{~}{r}+10(\stackrel{~}{L}a\stackrel{~}{E})^2\right)_{\mathrm{isco}}`$ (29)
$`=`$ $`{\displaystyle \frac{1}{1296}}\text{for }a=0,`$ (30)
$`\beta `$ $`=`$ $`{\displaystyle \frac{2}{\stackrel{~}{r}_{\mathrm{isco}}^4}}\left((\stackrel{~}{L}a^2\stackrel{~}{E}\stackrel{~}{\mathrm{\Omega }})\stackrel{~}{r}3(\stackrel{~}{L}a\stackrel{~}{E})(1a\stackrel{~}{\mathrm{\Omega }})\right)_{\mathrm{isco}}`$ (31)
$`=`$ $`{\displaystyle \frac{1}{36\sqrt{3}}}\text{for }a=0.`$ (32)
Here $`\stackrel{~}{r}_{\mathrm{isco}}`$ and $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}`$ are given by Eqs. (7) and (8); and $`\stackrel{~}{L}_{\mathrm{isco}}`$ and $`\stackrel{~}{E}_{\mathrm{isco}}`$ are expressed in terms of $`\stackrel{~}{r}_{\mathrm{isco}}`$ by Eqs. (9) and (10). Numerical values of $`\alpha `$ and $`\beta `$, computed from these equations, are tabulated in Table I.

### C Solution for Motion in the Transition Regime

The equation of motion in the transition regime, Eq. (26), can be converted into dimensionless form by setting

$$R=\eta ^{2/5}R_oX,\stackrel{~}{\tau }=\eta ^{1/5}\tau _oT,$$
(33)
where

$$R_o=(\beta \kappa )^{2/5}\alpha ^{3/5},\tau _o=(\alpha \beta \kappa )^{1/5};$$
(34)
cf. Table I. The resulting dimensionless equation of motion is

$$\frac{d^2X}{dT^2}=X^2T.$$
(35)
We seek the unique solution of this differential equation which, at early times $`T1`$, joins smoothly onto the adiabatic inspiral solution of Sec. II. In that adiabatic inspiral, the orbit is the circle at the minimum of the effective potential of Fig. 1 and Eq. (23), i.e., the circle at $`R=\sqrt{\beta \xi /\alpha }=\sqrt{\beta \kappa \eta \tau /\alpha }`$, which translates into

$$X=\sqrt{T}\text{for adiabatic inspiral near the isco.}$$
(36)
We have not been able to find an analytic formula for the solution to the equation of motion (35) that joins smoothly onto this adiabatic solution, but it is easy to construct the unique solution numerically. It is plotted in Figures 2 and 3, along with the adiabatic inspiral solution (36) and the plunge solution \[Eq. (38) below\].

The transition solution is well approximated by adiabatic inspiral at times $`T<1`$, but at $`T>1`$ it deviates from adiabatic inspiral and evolves smoothly into a plunge. The solution diverges ($`X\mathrm{}`$) at a finite time $`T=T_{\mathrm{plunge}}3.412`$.The divergence of $`X`$ at $`T=T_{\mathrm{plunge}}`$ does not imply a divergence of $`r`$ or any other physical quantity. Rather, it marks the breakdown of the transition approximation at very large values of $`|X|`$, $`|X|X_{\mathrm{break}}\eta ^{2/5}`$; \[cf. Eq. (33)\]. More specifically, when $`XX_{\mathrm{break}}`$, higher-order terms in the Taylor expansion (23) become important and stop the divergence. Well before this (in fact, throughout the range $`1XX_{\mathrm{break}}`$) both the transition approximation and the free-fall approximation ($`\stackrel{~}{E}=\stackrel{~}{E}_{\mathrm{final}}=\text{constant}`$, $`\stackrel{~}{L}=\stackrel{~}{L}_{\mathrm{final}}=\text{constant}`$) are valid, so these two approximations can be matched in this regime to obtain a solution valid all the way down to the horizon. The same type of breakdown also occurs at the other asymptotic limit $`+XX_{\mathrm{break}}`$: The transition regime’s adiabatic-inspiral equation (36) breaks down and must be replaced, via matching at $`1XX_{\mathrm{break}}`$, by the exact Kerr metric’s adiabatic inspiral formulae .

In the plunge regime, radiation reaction is unimportant; i.e., the orbit evolves inward with (very nearly) constant orbital angular momentum $`\stackrel{~}{L}\stackrel{~}{L}_{\mathrm{final}}`$ and energy $`\stackrel{~}{E}\stackrel{~}{E}_{\mathrm{final}}`$ (which we evaluate below); i.e., the orbit is well approximated by geodesic free fall. In the dimensionless equation of motion (35), the free-fall approximation translates into neglecting the last term, $`T`$, so $`d^2X/dT^2=X^2`$, which has the analytic first integral

$$dX/dT=\sqrt{\mathrm{constant}\frac{2}{3}X^3}.$$
(37)
For large $`|X|`$, the constant can be neglected and we obtain the analytic solution

$$X=\frac{6}{(T_{\mathrm{plunge}}T)^2}\text{for plunge near the isco,}$$
(38)
which is plotted in Figs. 2 and 3.

Combining Eqs. (16), (25), and (33), one finds that throughout the transition regime, the energy and angular-momentum deficits (i.e., the deviations of $`\stackrel{~}{E}`$ and $`\stackrel{~}{L}`$ from their isco values) scale as $`\eta ^{4/5}`$. In particular, the final deficits in the plunge stage are given by

$`\stackrel{~}{L}_{\mathrm{final}}\stackrel{~}{L}_{\mathrm{isco}}`$ $`=`$ $`(\kappa \tau _0T_{\mathrm{plunge}})\eta ^{4/5},`$ (39)
$`\stackrel{~}{E}_{\mathrm{final}}\stackrel{~}{E}_{\mathrm{isco}}`$ $`=`$ $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}(\kappa \tau _0T_{\mathrm{plunge}})\eta ^{4/5},`$ (40)
where, as was noted above,

$$T_{\mathrm{plunge}}=3.412.$$
(41)
## IV Gravitational Waves from Transition Regime, and their Observability

The gravitational waves emitted in the transition regime are all near the orbital frequency $`2\pi \mathrm{\Omega }_{\mathrm{isco}}`$ and its harmonics. The strongest waves are at the second harmonic (twice the orbital frequency):

$$f2\frac{\mathrm{\Omega }_{\mathrm{isco}}}{2\pi }=\frac{\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}}{\pi M}.$$
(42)
We shall compute their properties.

The transition waves last for a proper time $`\mathrm{\Delta }\tau =M\mathrm{\Delta }\stackrel{~}{\tau }=M\eta ^{1/5}\stackrel{~}{\tau }_o\mathrm{\Delta }T`$, during which the body spirals inward through a radial distance $`\mathrm{\Delta }r=M\mathrm{\Delta }R=M\eta ^{2/5}R_o\mathrm{\Delta }X`$, where $`\mathrm{\Delta }T`$ covers the range $`T1`$ to $`2.3`$ and $`\mathrm{\Delta }X`$ covers the range $`X1`$ to $`X5`$ (Fig. 3); i.e.,

$$\mathrm{\Delta }T=3.3,\mathrm{\Delta }X=6.$$
(43)
Correspondingly, neglecting any cosmological redshift, the duration of the transition waves as seen at Earth is

$$\mathrm{\Delta }t=\frac{M}{(d\stackrel{~}{\tau }/d\stackrel{~}{t})_{\mathrm{isco}}}\eta ^{1/5}\stackrel{~}{\tau }_o\mathrm{\Delta }T,$$
(44)
and their frequency band is $`\mathrm{\Delta }f=(1/\pi M)(d\stackrel{~}{\mathrm{\Omega }}/d\stackrel{~}{r})_{\mathrm{isco}}\mathrm{\Delta }\stackrel{~}{r}`$, which, using the above expression for $`\mathrm{\Delta }r`$ and Eq. (1) for $`\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{r})`$, gives

$$\mathrm{\Delta }f=\frac{3}{2\pi M}\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}^2\sqrt{\stackrel{~}{r}_{\mathrm{isco}}}\eta ^{2/5}R_o\mathrm{\Delta }X.$$
(45)
The total number of cycles of these transition waves is

$$N_{\mathrm{cyc}}=f\mathrm{\Delta }t=\frac{\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}\stackrel{~}{\tau }_o}{\pi (d\stackrel{~}{\tau }/d\stackrel{~}{t})_{\mathrm{isco}}}\eta ^{1/5}\mathrm{\Delta }T.$$
(46)
These second-harmonic waves arriving at Earth have the form $`h_+=h_{+\mathrm{amp}}\mathrm{cos}(2\pi f𝑑t+\phi _+)`$, $`h_\times =h_{\times \mathrm{amp}}\mathrm{cos}(2\pi f𝑑t+\phi _\times )`$, where $`\phi _+`$ and $`\phi _\times `$ are constant phases. The amplitudes $`h_{+\mathrm{amp}}`$ and $`h_{\times \mathrm{amp}}`$ depend on the source’s orientation. When one squares and adds these amplitudes and then averages over the sky (“$`\mathrm{}`$”), one obtains an rms amplitude

$$h_{\mathrm{amp}}^{\mathrm{rms}}=h_{+\mathrm{amp}}^2+h_{\times \mathrm{amp}}^2^{1/2},$$
(47)
which is related to the power being radiated into the second harmonic by $`\dot{E}_2=4\pi D^2(h_{\mathrm{amp}}^{\mathrm{rms}})^2(2\pi f)^2/(32\pi )`$; cf. Eq. (35.27) of MTW . Here $`D`$ is the distance to the source. Equating this to the radiated power $`\dot{E}_2=(32/5)\eta ^2\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}^{10/3}\dot{}_{\mathrm{},2}`$ , where $`\dot{}_{\mathrm{},2}`$ is a relativistic correction factor listed on the first line of Table IV of , we obtain the following expression for the waves’ rms amplitude

$$h_{\mathrm{amp}}^{\mathrm{rms}}=\frac{8}{\sqrt{5}}\frac{M\eta }{D}\stackrel{~}{\mathrm{\Omega }}_{\mathrm{isco}}^{2/3}\sqrt{\dot{}_{\mathrm{},2}}.$$
(48)
The signal to noise ratio $`S/N`$ that these waves produce in LISA depends on the orientations of LISA and the source relative to the line of sight between them. When one squares $`S/N`$ and averages over both orientations, then takes the square root, one obtains

$$\left(\frac{S}{N}\right)_{\mathrm{rms}}=\frac{h_{\mathrm{amp}}^{\mathrm{rms}}}{\sqrt{5S_h(f)/\mathrm{\Delta }t}}.$$
(49)
Here $`5S_h(f)`$ is the spectral density of LISA’s strain noise inverse-averaged over the sky§§§ i.e., $`1/(5S_h)`$ average over the sky of 1/(spectral density). The factor $`5`$ in this definition is to produce accord with the conventional notation for ground-based interferometers, where $`S_h(f)`$ denotes the spectral density for waves with optimal direction and polarization. In the case of LISA, at frequencies above about 0.01 Hz, the beam pattern shows sharp frequency-dependent variations with direction due to the fact that the interferometer arms are acting as one-pass delay lines rather than optical cavities, and this produces a more complicated dependence of sensitivity on angle than for ground-based interferometers. As a result, $`S_h`$ (as we have defined it) is the spectral density for optimal direction and polarization only below about 0.01 Hz, not above.  and $`1/\mathrm{\Delta }t`$ is the band width associated with the waves’ duration $`\mathrm{\Delta }t`$.

The noise spectral density $`S_h(f)`$ for the current straw-man design of LISA has been computed by the LISA Mission Definition Team . An analytic fit to this $`S_h(f)`$, after averaging over some small-amplitude oscillations that occur at $`f>0.01`$ Hz, is the following:

$`S_h(f)=`$ $`[(4.6\times 10^{21})^2+(3.5\times 10^{26})^2\left({\displaystyle \frac{1\mathrm{H}\mathrm{z}}{f}}\right)^4`$ (51)
$`+(3.5\times 10^{19})^2\left({\displaystyle \frac{f}{1\mathrm{H}\mathrm{z}}}\right)^2]\mathrm{Hz}^1.`$
The rate for $`\mu 10M_{}`$ black holes to spiral into $`M10^6M_{}`$ black holes in galactic nuclei has been estimated by Sigurdsson and Rees ; their “very conservative” result is $``$ one event per year out to 1 Gpc. The inspiraling holes are likely to be in rather eccentric, nonequatorial orbits , for which our analysis needs to be generalized. If, however, the orbit is circular and equatorial and the holes are at 1 Gpc distance, then the above formulas give the numbers shown in Table II.

As shown in the table, the signal to noise for this source is of order unity. With some luck in the orientation of LISA, the orientation of the source, the distance to the source, and/or the holes’ masses, a $`S/N`$ of a few might occur. Since the signal would already have been detected from the much stronger adiabatic inspiral waves, this signal strength could be enough to begin to explore the details of the transition from inspiral to plunge.

## V Conclusions

Our analysis of the transition regime has been confined to circular, equatorial orbits. This is a serious constraint, since there is strong reason to expect that most inspiraling bodies will be in orbits that are strongly noncircular and nonequatorial . Our estimated signal-to-noise ratio, $`S/N1`$, for LISA’s observations of the transition regime from a circular, equatorial orbit at the plausible distance $`1`$ Gpc suggests that for more realistic orbits the transition regime might be observable. This prospect makes it important to generalize our analysis to more realistic orbits.

Full analyses for equatorial, noncircular orbits and for nonequatorial, circular orbits can be carried out using techniques now in hand: the Teukolsky formalism, and computations of the orbital evolution based on the energy and angular momentum radiated down the hole and off to infinity (see, e.g., Ref.  and references therein). For nonequatorial, noncircular orbits, the analysis should also be possible with existing techniques — up to an unknown radiation-reaction-induced rate of evolution of the Carter constant. That unknown quantity could be left as a parameter in the analysis, to be determined when current research on gravitational radiation reaction  has reached fruition.

When this paper was in near final form, we became aware of a similar analysis, by Buonanno and Damour , of the transition from inspiral to plunge. Whereas we treat the case of infinitesimal mass ratio $`\eta 1`$ and finite black-hole spin $`1<a<+1`$, Buananno and Damour treat finite $`\eta `$ ($`0<\eta 1/4`$) and vanishing spins $`a=0`$. Both analyses give the same dimensionless equation of motion (35) for the transition regime.

## Acknowledgments

We thank James Anderson for helpful discussions, Theocharis Apostolatos and Richard O’Shaughnessy for checking some details of our analysis, and Sam Finn for helpful advice and for permission to use unpublished results from his numerical solutions of the Teukolsky equation; see Tables I and II. This work was supported in part by NASA grant NAG5-6840, and in view of its potential applications to LIGO, also by NSF grant AST-9731698.


format: text

text:
# UCLA/00/TEP/5 MPI-PhT/2000-05 March 2000 Measurement of the gluon PDF at small 𝑥 with neutrino telescopes

## I Introduction

Atmospheric neutrinos and muons are the most important source of background for present and future neutrino telescopes, which are expected to open a new window in astronomy by detecting neutrinos from astrophysical sources.

At energies above 1 TeV, atmospheric lepton fluxes have a prompt component consisting of neutrinos and muons created in semileptonic decays of charmed particles, as opposed to the conventional leptons coming from decays of pions and kaons. Thus a model for charm production and decays in the atmosphere is required.

We base our model on QCD, the theoretically preferred model, to compute the charm production. We use a next-to-leading order perturbative QCD (NLO pQCD) calculation of charm production in the atmosphere, followed by a full simulation of particle cascades generated with PYTHIA routines .

We have already examined the prompt muon and neutrino fluxes in two previous papers  (called GGV1 and GGV2 from now on).

In our first paper , we found that the NLO pQCD approach produces fluxes in the bulk of older predictions (not based on pQCD) as well as of the recent pQCD semianalytical analysis of Pasquali, Reno and Sarcevic . We also explained the reason of the low fluxes of the TIG model , the first to use pQCD in this context, which were due to the chosen extrapolation of the gluon partonic distribution function (PDF) at small momentum fractions $`x`$, and we confirmed the overall validity of their pioneering approach to the problem.

In our second paper , we analyzed in detail the dependence of the fluxes on the extrapolation of the gluon PDF at small $`x`$, which, according to theoretical models, is assumed to be a power law with exponent $`\lambda `$,

$$xg(x)x^\lambda ,$$
(1)
with $`\lambda `$ in the range 0–0.5. Particle physics experiments are yet unable to determine the value of $`\lambda `$ at $`x<10^5`$. We found that the choice of different values of $`\lambda `$ at $`x<10^5`$ leads to a wide range of final background fluxes at energies above 105 GeV.

Due to this result, in GGV2 we suggested the possibility of measuring $`\lambda `$ through the atmospheric leptonic fluxes at energies above $`10^5`$ GeV, not the absolute fluxes, because of their large theoretical error, but rather through their spectral index (i.e. the “slope”). In particular, we now propose to use the slope of the flux of down-going muons.

We want to stress that we are proposing to use down-going muons, at energies $`E_\mu 100`$ TeV, where prompt muons dominate over conventional ones, and not up-going neutrino-induced muons whose flux is orders of magnitude smaller. While an important contribution to up-going muons is expected from astrophysical neutrinos, there is no background for down-going atmospheric muons.

In this paper we further investigate the possibility of measuring $`\lambda `$, in the more general context of an overall error analysis of our model.

We can identify five potential causes of uncertainty in our final results. The first one is the possible presence of large logarithms of the type $`\alpha _s\mathrm{ln}p_T^2`$ and $`\alpha _s\mathrm{ln}s`$ (the latter are the so called “$`\mathrm{ln}(1/x)`$” terms). The second is the treatment of the multiplicity in the production of $`c\overline{c}`$ at high energies.

The third one consists of all the sources of uncertainty hidden in the treatment of particle cascades generated by PYTHIA. The fourth one is the uncertainty in the NLO pQCD charm production model we use. This includes the dependence of the fluxes on the three parameters of the model and the PDF’s used. The fifth and final one is the choice of the primary cosmic ray flux, which is the input of our simulation. Of all these potential sources of uncertainty we conclude that only the last two are relevant.

We deal with these five potential sources of error in turn. In Sect. II we address the question of the large logarithms $`\alpha _s\mathrm{ln}p_T^2`$ and $`\alpha _s\mathrm{ln}s`$, and in Sect. III we analyze the problem of multiplicity in our charm production mechanism.

In Sect. IV we consider the uncertainties due to the cascade generation by PYTHIA and to our NLO pQCD charm production (the core of our analysis). Here we evaluate the errors due to the parameters of the model, errors that affect the charm production cross section, the final differential and integrated fluxes and their spectral indices. We also determine how the final results (fluxes and their spectral indices) are affected by the choice of different extrapolations of the PDF’s at $`x<10^5`$.

We are finally ready in Sect. V to discuss how $`\lambda `$ could be measured. We study the dependence on the different extrapolations of $`\lambda `$ at $`x<10^5`$, we consider the spectral indices and, using the discussion of Sect. IV, we provide an estimate of the errors on these indices and examine the feasibility of an experimental determination of $`\lambda `$ at $`x<10^5`$ with neutrino telescopes.

Finally in Sect. VI we discuss the error on the determination of $`\lambda `$ coming from the uncertainties in the elemental composition of the cosmic ray flux.

The determination of $`\lambda `$ at small $`x<10^5`$ is important because in this range saturation, unitarity and shadowing effects should become important. The PDF sets we use have been derived without including saturation effects. Even if this procedure seems to work very well in the HERA regime (where there might be some indications of saturation already, see e.g. ), here we are extrapolating the resulting gluon PDF’s at even smaller $`x`$ values where saturation may become important. With respect to unitarity, using the expression of the Froissart bound on the gluon structure functions given in Eq. 31 of Ref.  we see that the extrapolated gluon PDF’s we use, with $`\lambda =0.40.5`$, violate this bound at $`x`$ values between 0.5 $`\times 10^7`$ and 1 $`\times 10^7`$, for the characteristic momenta $`Q^2m_c^23`$ GeV2 we have, which corresponds to leptonic energies of $`12\times 10^6`$ GeV. Always using the Froissart bound on the gluon PDF as given in the Eq. 31 of Ref. , the gluon PDF‘s extrapolated with $`\lambda 0.3`$ encounter this bound at $`x<10^8`$, what corresponds to leptonic energies larger than $`10^8`$ GeV, beyond the energy range relevant in this paper. Shadowing of the gluons in the atmospheric nucleons and nuclei, which we have not included here, could decrease the amplitude of the gluon PDF’s by about $`10\%`$ at $`x10^5`$ and up to as much as $`30`$ to $`40\%`$ at $`x10^7`$ (see e.g. ), what would also change the effective value of $`\lambda `$. There are no shadowing effects in the cosmic ray nucleons and nuclei. Only the dominant $`x`$ of the gluons in the atmosphere is small in our calculation, while the dominant $`x`$ of the partons in the cosmic rays is large. Thus shadowing effects do not depend on the unknown composition of cosmic rays, but could only be important for the atmospheric partons. The reason for the different characteristic values of $`x`$ in the target and projectile partons is the following (for more details see GGV2). Due to the steep decrease with increasing energy of the incoming flux of cosmic rays, only the most energetic charm quarks produced count and those come from the interaction of projectile partons carrying a large fraction of the incoming nucleon momentum. Thus, the characteristic $`x`$ of the projectile parton, $`x_1`$, is large. It is $`x_1O(10^1)`$. We can, then, inmediately understand that very small parton momentum fractions are needed in our calculation, because typical partonic center of mass energies $`\sqrt{\widehat{s}}`$ are close to the $`c\overline{c}`$ threshold, $`2m_c2`$ GeV, (since the differential $`c\overline{c}`$ production cross section decreases with increasing $`\widehat{s}`$) while the total center of mass energy squared is $`s=2m_NE`$ (with $`m_N`$ the nucleon mass, $`m_N1`$ GeV). Calling $`x_2`$ the momentum fraction of the target parton (in a nucleous of the atmosphere), then, $`x_1x_2\widehat{s}/s=4m_c^2/(2m_NE)`$ GeV/$`E`$. Thus, $`x_2O`$(GeV/0.1 $`E`$), where $`E`$ is the energy per nucleon of the incoming cosmic ray in the lab. frame. The characteristic energy $`E_c`$ of the charm quark and the dominant leptonic energy $`E_l`$ in the fluxes are $`E_lE_c0.1E`$, thus $`x_2O`$(GeV/ $`E_l`$).

## II Importance of the $`\alpha _s\mathrm{ln}(1/x)`$ terms

We address here a concern that has been expressed to us several times, about the applicability of perturbative QCD calculations, mostly done for accelerator physics, to the different kinematic domain of cosmic rays.

Contrary to the case in accelerators, we do not have the uncertainty present in the differential cross sections  when the transverse momentum $`p_T`$ is much larger than $`m_c`$, uncertainty which is due to the presence of large logarithms of $`(p_T^2+m_c^2)/m_c^2`$. The reason is that we do not have a forward cut in acceptance, and so the characteristic transverse charm momentum in our simulations is of the order of the charm mass, $`p_TO(m_c)`$, and not $`p_TO(m_c)`$ as in accelerator experiments.

We may however, depending on the steepness of the gluon structure function $`\lambda `$, have large logarithms of the type $`\alpha _s\mathrm{ln}s`$, known as “$`\mathrm{ln}(1/x)`$” terms (here $`x\sqrt{4m_c^2/s}`$ is the average value of the hadron energy fraction needed to produce the $`c\overline{c}`$ pair at hadronic center of mass energy squared $`s`$). These “$`\mathrm{ln}(1/x)`$” terms arise when the t-channel gluon exchange becomes large, and eventually they have to be resummed. Although techniques exist for resumming these logarithms , we have not done it. On the other hand we have phenomenologically altered the behavior of the parton distribution functions at small $`x`$ by imposing a power law dependence of the form $`xf(x)x^\lambda `$. This is analogous to resumming the $`\mathrm{ln}(1/x)`$ terms in a universal fashion and absorbing them in an improved evolution equation for the gluon density (such as the Balitskyiĭ-Fadin-Kuraev-Lipatov (BFKL) evolution equation) , a procedure which increases $`\lambda `$. For sufficiently large $`\lambda `$, the large $`\mathrm{ln}(1/x)`$ terms should not be present.

To find if our NLO $`c\overline{c}`$ cross sections are dominated by the $`\mathrm{ln}(1/x)`$ terms, we have used the following qualitative criterion . We have plotted the ratio

$$R=\frac{\sigma _{NLO}\sigma _{LO}}{\sigma _{LO}\alpha _s\mathrm{ln}(s/m_c^2)/\pi }$$
(2)
as a function of the beam energy $`E`$. If the ratio is constant we are dominated by the $`\mathrm{ln}(1/x)`$ terms and if it decreases we are not. The good behavior is a decreasing $`R`$. Figure 1 shows indeed that up to highest energy we consider in this paper, i.e. $`10^{11}`$ GeV, $`R`$ decreases for $`\lambda 0.2`$, but is roughly constant for smaller $`\lambda `$’s. This indicates that we are not dominated by the $`\mathrm{ln}(1/x)`$ terms provided $`\lambda 0.2`$.

Clearly, this test involving the $`R`$ ratio does not say anything about $`ln(1/x)`$ higher order corrections. One can only argue that if the $`ln(1/x)`$ terms are not dominant at NLO (for $`R`$ decreasing with energy) the corresponding $`[ln(1/x)]^n`$ terms may also be non-dominant in higher order corrections. In any event, the data on charm production that could be inferred at $`x<10^5`$, from the slope of atmospheric muon fluxes, really give information on the product of the gluon PDF and the parton cross section and a measurement of this product is useful. One can expect that the $`ln(1/x)`$ terms at higher order may be better understood by the time the data will come.

## III Multiplicity in charm production

Another concern is the fact that at high energies the charm production cross section we use is sometimes larger than the total $`pp`$ cross section. At first sight this seems absurd, but we show here that the cross section we use is the inclusive cross section, which contains the charm multiplicity, i.e. it counts the number of $`c\overline{c}`$ pairs produced, and so can be larger than the total cross section. On the other hand, the contribution of $`c\overline{c}`$ producing events to the total $`pp`$ cross section, i.e. the cross section for producing at least one $`c\overline{c}`$ pair, is always smaller than the total $`pp`$ cross section.

We call $`\sigma _{QCD}`$ the perturbative QCD cross section of $`c\overline{c}`$ pair production in $`pp`$ collisions,

$$\sigma _{QCD}=\underset{ij}{}\sigma _{QCD}(ijc\overline{c}),$$
(3)
where the sum is over the partons $`i`$ and $`j`$ in the colliding nucleons, and

$$\sigma _{QCD}(ijc\overline{c})=𝑑x_1𝑑x_2𝑑Q^2\frac{d\widehat{\sigma }(ijc\overline{c})}{dQ^2}f_i(x_1,\mu _F^2)f_j(x_2,\mu _F^2).$$
(4)
Here $`d\widehat{\sigma }(ijc\overline{c})/dQ^2`$ is the $`ijc\overline{c}`$ parton scattering cross section, $`Q^2`$ is the four-momentum transfer squared, $`x_i`$ is the fraction of the momentum of the parent nucleon carried by parton $`i`$, and $`f_i(x,\mu _F^2)`$ is the usual parton distribution function for parton momentum fraction $`x`$ and factorization scale $`\mu _F`$.

In the scattering of each pair of partons (one parton from the target and one from the projectile) only one $`c\overline{c}`$ pair may be produced, but the number of parton pairs interacting in each nucleon-nucleon collision is in general not limited to one and it increases with the number of partons $`f(x,\mu _F^2)dx`$ in each nucleon.

For $`\lambda `$ close to 0.5, $`\sigma _{QCD}`$ becomes larger than the total $`pp`$ cross section $`\sigma _{pp}200\mathrm{mb}`$ at $`E_p10^{10}`$ GeV. It is obvious therefore that our results at high energy and large $`\lambda `$ are unphysical, unless multiplicity is taken into account. In fact, multiparton interactions should be taken into account already at a smaller cross section of order 10 mb, as determined in studies of double parton scattering .

In order to incorporate multiparton scatterings into our analysis, we use an impact-parameter representation for the scattering amplitude, and ignore spin-dependent effects (cfr. ). Assuming the validity of factorization theorems, the mean number of parton-parton interactions $`ijc\overline{c}`$ in the collision of two protons at impact parameter $`\stackrel{}{b}`$ is given by

$$n_{c\overline{c}}(\stackrel{}{b})=\underset{ij}{}d^2b^{}𝑑x_1𝑑x_2𝑑Q^2\frac{d\widehat{\sigma }(ijc\overline{c})}{dQ^2}f_i(x_1,\mu _F^2,\stackrel{}{b^{}})f_j(x_2,\mu _F^2,\stackrel{}{b}+\stackrel{}{b^{}}),$$
(5)
where $`f_i(x,\mu _F^2,\stackrel{}{b})d^2b`$ is the number of partons $`i`$ in the interval ($`x`$, $`x+dx`$) and in the transverse area element $`d^2b`$ at a distance $`\stackrel{}{b}`$ from the center of the proton. For simplicity of notation we drop the vector symbol in $`\stackrel{}{b}`$ and write $`b`$ from now on.

If $`n_{c\overline{c}}(b)1`$, $`n_{c\overline{c}}(b)`$ is the probability of producing a $`c\overline{c}`$ pair in a $`pp`$ collision at impact parameter $`b`$. If $`n_{c\overline{c}}(b)1`$, $`n_{c\overline{c}}(b)`$ is just the mean value of $`k`$, the number of $`c\overline{c}`$ pairs produced, at impact parameter $`b`$. Let the probability of $`k`$ scatterings $`ijc\overline{c}`$ in a $`pp`$ collision at impact parameter $`b`$ be $`P_{kc\overline{c}}(b)`$. Then

$$n_{c\overline{c}}(b)=\underset{k=0}{\overset{\mathrm{}}{}}kP_{kc\overline{c}}(b).$$
(6)
The $`k`$-tuple parton cross section is obtained by integrating the probability of exactly $`k`$ parton scatterings $`P_{kc\overline{c}}(b)`$ over the impact parameter $`b`$,

$$\sigma _{kc\overline{c}}=d^2bP_{kc\overline{c}}(b),$$
(7)
the inclusive cross section for charm production is, thus, $`\sigma _{c\overline{c}\mathrm{incl}}=_kk\sigma _{kc\overline{c}}`$ and the contribution of charm producing processes to the total cross section is $`\sigma _{c\overline{c}}=_k\sigma _{kc\overline{c}}`$ for $`k0`$.

In our evaluation of charm production by cosmic ray interactions in the atmosphere, we must count the number of $`c\overline{c}`$ pairs produced in the $`pp`$ collision. So we are precisely interested in the inclusive cross section $`\sigma _{c\overline{c}\mathrm{incl}}`$, which includes the number $`k`$ of $`c\overline{c}`$ pairs produced per collision (the multiplicity). We find

$$\sigma _{c\overline{c}\mathrm{incl}}=\underset{k}{}k\sigma _{kc\overline{c}}=d^2b\underset{k}{}kP_{kc\overline{c}}(b)=d^2bn_{c\overline{c}}(b).$$
(8)
This cross section can be larger than the total $`pp`$ cross section, because it accounts for multiparton interactions. In particular, using $`\sigma _{c\overline{c}}`$, the contribution of charm producing processes to the total cross section defined above, the ratio $`\sigma _{c\overline{c}\mathrm{incl}}/\sigma _{c\overline{c}}`$ gives the average charm multiplicity.

Notice that here we consider only independent production of $`c\overline{c}`$ pairs, so that from each pair of colliding partons it results only one $`c\overline{c}`$ pair, and we neglect coherent production of multiple $`c\overline{c}`$ pairs in 2$``$4, 2$``$6, etc. processes. This will underestimate the charm production cross section.

We assume in the following that the partonic distributions $`f_i(x,\mu _F^2,b)`$ factorize as

$$f_i(x,\mu _F^2,b)=f_i(x,\mu _F^2)\rho _i(b),$$
(9)
where $`f_i(x,\mu _F^2)`$ is the usual parton distribution function, and $`\rho _i(b)`$ is the probability density of finding a parton in the area $`d^2b`$ at impact parameter $`b`$. We normalize $`\rho _i(b)`$ to $`d^2b\rho _i(b)=1`$, to maintain the usual normalization $`𝑑xxf_i(x)=1`$. The factorization in Eq. (9) is consistent with the usual parton picture and with our assumption of no parton-parton correlations.

The mean number of $`ijc\overline{c}`$ scatterings at impact parameter $`b`$ then becomes

$$n_{c\overline{c}}(b)=\underset{ij}{}a_{ij}(b)\sigma _{QCD}(ijc\overline{c}),$$
(10)
where

$$a_{ij}(b)=d^2b^{}\rho _i(b^{})\rho _j(b+b^{})$$
(11)
is an overlap integral, and $`\sigma _{QCD}(ijc\overline{c})`$ is the QCD parton-parton cross section for $`ijc\overline{c}`$, as in Eq. (4). From the normalization of $`\rho _i(b)`$ it follows that $`d^2ba_{ij}(b)=1`$ for every $`i,j`$. Hence from Eqs. (8) and (10) we find

$$\sigma _{c\overline{c}\mathrm{incl}}=\sigma _{QCD},$$
(12)
where $`\sigma _{QCD}`$, given is Eq. (3), is the charm production cross section calculated within QCD. This justifies our use of $`\sigma _{QCD}`$ as $`\sigma _{c\overline{c}\mathrm{incl}}`$ in the calculation of the atmospheric fluxes.

The way in which we use $`\sigma _{c\overline{c}\mathrm{incl}}`$ in our simulation is as follows. We input only one $`c\overline{c}`$ pair per $`pp`$ collision at a given energy $`E`$, and multiply the outcome by $`\sigma _{c\overline{c}\mathrm{incl}}`$, which includes the $`c\overline{c}`$ multiplicity. We make, therefore, the following approximation in the kinematics of the $`c\overline{c}`$ pairs produced in the same $`pp`$ interaction. Even if in a real multiparton collision the energy available to the second and other $`c\overline{c}`$ pairs is smaller than $`E`$, we are neglecting this difference. This is a very good approximation because the fraction of center of mass energy that goes into a $`c\overline{c}`$ pair is of the order of $`\sqrt{\widehat{s}/s}\sqrt{10\mathrm{GeV}/E}1`$ at the high energies we are concerned with.

Although we have explicitly proven Eq. (12) in the absence of parton–parton correlations, the same result can be obtained when correlations are present (see sect. V of Ref.  and references therein). What is proven even in the presence of correlations is that the pQCD single scattering cross section $`\sigma _{QCD}`$ is equal to the average number of parton-parton collisions, call it $`<N>`$, multiplied by the contribution of $`c\overline{c}`$ producing events to the total cross section (the cross-section $`\sigma _{c\overline{c}}`$ defined above), namely $`\sigma _{QCD}=<N>\sigma _{c\overline{c}}`$ (while the inclusive cross section is equal to the average multiplicity of $`c\overline{c}`$ pairs multiplied by $`\sigma _{c\overline{c}}`$). $`<N>`$ may in general contain contributions from two types of collisions. One type consists of collisions of pairs of partons (consisting of one parton from each interacting nucleon) which interact only once at different points of the transverse plane. Each collision of this type results in our case in one $`c\overline{c}`$ pair produced. The second type consists of rescatterings in which one parton of one of the nucleons and its interaction products interact with several partons of the other nucleon. In interactions of this type, which are much rarer than the first ones, the number of $`c\overline{c}`$ pairs produced not necessarily coincides with the number of collisions. If rescatterings can be neglected, then $`<N>`$ is the average number of $`c\overline{c}`$ pairs produced and $`\sigma _{QCD}`$ is the inclusive $`c\overline{c}`$ production cross section as stated in Eq. (12). Otherwise small rescattering corrections, to our knowledge not yet calculated , are necessary (rescatterings would also modify the energy spectrum of the particles produced).

## IV Uncertainties due to cascade simulation, parameters of charm production model and choice of PDF’s

In our first paper (GGV1) we considered the uncertainties related to the cascade generation in PYTHIA. There we tried different modes of cascade generation, different options allowed by PYTHIA in the various stages of parton showering, hadronization, interactions and decays, etc., without noticing substantial changes in the final results (differing at most by $`10\%`$). These uncertainties are however very difficult to quantify, due to the nature of the PYTHIA routines. Since these uncertainties are small, we neglect them in this analysis and continue to use PYTHIA with the options described in GGV1 as our main choice for the simulation: ‘single’ mode with showering, ‘independent’ fragmentation, interactions and semileptonic decays according to TIG.

In our ‘single’ mode we enter only one $`c`$ quark in the particle list of PYTHIA, and we multiply the result by a factor of $`2`$ to account for the initial $`\overline{c}`$ quark. PYTHIA performs the showering, standard independent fragmentation, and follows all the interactions and decays using default parameters and options. In GGV1 we have argued that this ‘single’ approach is equivalent to what we called ‘double’ mode, in which both $`c\overline{c}`$ partons are placed in the initial event list, in the first step of a standard cascade evolution. The ‘single’ option is chosen thus, because it reduces considerably the computing time.

Important sources of uncertainty are contained in our charm production model, which is NLO pQCD as implemented in the MNR program , calibrated at low energies.

The calibration procedure consisted in the following:

$``$ choosing a PDF set from those available and fixing the related value of $`\mathrm{\Lambda }_{QCD}`$;***We note that $`\mathrm{\Lambda }_{QCD}`$ can be chosen in the MNR program independently of the PDF and therefore can constitute an additional independent parameter of our model. We have opted however to choose the value of $`\mathrm{\Lambda }_{QCD}`$ assumed in the PDF set being used, for consistency.

$``$ choosing $`m_c`$, $`\mu _F`$ and $`\mu _R`$, which are the charm mass, the factorization scale and the renormalization scale respectively, so as to fit simultaneously both the total and differential cross sections from existing fixed target charm production experiments  at the energy of 250 GeV, without additional normalization factors;

$``$ checking that the total cross section generated after the previous choices fits reasonably well the other existing experimental points for fixed target charm production experiments .

Besides the choice of the PDF set, our procedure has the freedom to choose reasonable values of the three parameters $`m_c`$, $`\mu _F`$, and $`\mu _R`$ so as to fit the experimental data. In GGV1 and GGV2 we made the standard choice  of

$$\mu _F=2m_T,\mu _R=m_T,$$
(13)
where $`m_T=\sqrt{p_T^2+m_c^2}`$ is the transverse mass. The values of the charm mass are taken slightly different for each PDF set, namely:

$`m_c`$ $`=`$ $`1.185\mathrm{GeV}\text{for MRS R1,}`$ (14)
$`m_c`$ $`=`$ $`1.310\mathrm{GeV}\text{for MRS R2,}`$ (15)
$`m_c`$ $`=`$ $`1.270\mathrm{GeV}\text{for CTEQ 4M,}`$ (16)
$`m_c`$ $`=`$ $`1.250\mathrm{GeV}\text{for MRST.}`$ (17)
Here we explore the changes induced in cross sections and fluxes at high energies by different choices of $`m_c`$, $`\mu _F`$, and $`\mu _R`$ which fulfil our calibration requirements.

We have performed this analysis with the most recent PDF set: the MRST  (other PDF’s give similar results). At first we fix $`\lambda =0`$ and then we examine other values of $`\lambda `$. We note that the calibration procedure described above is independent of $`\lambda `$ because it involves only relatively low energies, where the low $`x`$ extrapolation is not an issue.

### A MRST $`\lambda =0`$: fluxes

We considered the $`\lambda =0`$ case because it is the most significant for the evaluation of the uncertainties in the spectral indices, as it will be clear in the next subsection. We have considered the following reasonable ranges of the parameters

$`1.1\mathrm{GeV}m_c1.4\mathrm{GeV},`$ (18)
$`0.5m_T\mu _F2m_T,`$ (19)
$`0.5m_T\mu _R2m_T,`$ (20)
where the bounds on $`m_c`$ come from the 1998 Review of Particle Physics , while those for $`\mu _F`$ and $`\mu _R`$ are those used in the existing literature .

Within these ranges we have looked for values of the three parameters capable of reproducing the experimental data in our calibration procedure described before. Table I summarizes the different sets of parameters: we have varied the charm mass through the values $`m_c=1.1`$, 1.2, 1.25, 1.3, 1.4 GeV ($`m_c=1.25`$ GeV was our previous optimal choice for MRST in Eq. (17)) and then, for each value of $`m_c`$, we have found values of the factorization and renormalization scales that reproduce the experimental total cross section $`\sigma _{c\overline{c}}=13.5\pm 2.2`$ $`\mu `$b at 250 GeV . In particular, for each value of $`m_c`$, we took $`\mu _F=m_T/2`$, $`m_T`$, $`2m_T`$ and then, for each $`m_c`$, $`\mu _F`$ pair, found the value of $`\mu _R`$ which best fits the data (see Table I).

We have also checked that these choices give good fits to the differential, besides the total, cross sections at 250 GeV , without additional normalization factors, as done for the original choice of parameters in GGV1. For $`m_c=1.1`$ GeV we had to choose values of $`\mu _R`$ slightly outside the range in Eq. (20) (but we have kept these values in our analysis anyway).

For all the sets of parameters in Table I we have run our full simulations for the MRST, $`\lambda =0`$ case and the results are described in Figs. 2-4.

In Figs. 2a and 2b we show the resulting total charm production cross section $`\sigma _{c\overline{c}}`$ for all of the fifteen sets of parameters in Table I, together with recent experimental data (from Table 1 of Ref. , where all the data for $`pp`$ or $`pN`$ collisions have been transformed into a $`\sigma _{c\overline{c}}`$ cross section following the procedure described in GGV1). Fig. 2b is an enlargement of the region of Fig. 2a containing the experimental data.

In Fig. 2a we see the spread of the cross sections, which is more than one order of magnitude at $`10^{11}`$ GeV. Above $`250GeV`$, one can clearly distinguish three “bands” of increasing cross sections for $`\mu _F=m_T/2`$, $`m_T`$ and $`2m_T`$. Within each “band” the cross sections increase with increasing values of $`m_c`$ (and correspondingly smaller values of $`\mu _R`$), in Table I. Our standard choice ($`m_c=1.25`$ GeV, $`\mu _F=2m_T`$, $`\mu _R=m_T`$) proves to be one of the highest cross sections we obtain.

In Fig. 2b we see better how all of these cross sections verify our calibration procedure. They pass through the point at 250 GeV , agree with the point at 400 GeV  and disagree with the very low experimental point at 200 GeV . The lower values of $`\mu _F=m_T/2`$ and $`m_T`$ fit better the lowest experimental point at 800 GeV , while the higher value of $`\mu _F=2m_T`$ fits better the upper point at 800 GeV .

We believe that the spread of the total cross sections shown in Fig. 2a provides a reasonable estimate of the uncertainty of our charm production model at fixed $`\lambda `$. Since our standard choice of parameters ($`m_c=1.25`$ GeV, $`\mu _F=2m_T`$, $`\mu _R=m_T`$) gives one of the highest cross sections (in better agreement with the more recent value of the cross section at 800 GeV ), the uncertainty band should be added under each of the cross section curves calculated with our standard choice of parameters (like the curves shown in Fig. 1 of GGV2).

Fig. 3 illustrates the corresponding spread of the final prompt fluxes. Although our results are for the MRST PDF’s extrapolated with $`\lambda =0`$ (the value of $`\lambda `$ which gives the lowest fluxes) similar spreads result from other PDF’s and $`\lambda `$’s. We show the flux of muons; the fluxes of muon-neutrinos and electron-neutrinos are essentially the same.

Similarly to what happens with cross sections in Fig. 2, the fluxes in Fig. 3 increase with $`\mu _F=m_T/2`$, $`m_T`$ and $`2m_T`$ and, within each band, they increase with increasing $`m_c`$ (and correspondingly smaller values of $`\mu _R`$), in Table I. At energies around $`10^6`$ GeV the total uncertainty is almost one order of magnitude and decreases slightly for higher energies. If we would decide to work only with $`\mu _F=2m_T`$ (which fits the experimental measurement at 800 GeV with the highest cross section), the uncertainty would be greatly reduced: the fluxes in this rather narrow band differ by less than 40%. We observe that the flux calculated with our standard choice of parameters ($`m_c=1.25`$ GeV, $`\mu _F=2m_T`$, $`\mu _R=m_T`$) is almost the highest, as was the case for the corresponding cross section in Fig. 2.

In Fig. 3 we also indicate the conventional and prompt fluxes from TIG; we notice that the TIG prompt flux is within our band of uncertainty, which is reasonable since TIG used a low $`\lambda =0.08`$ value for their predictions (see the discussions in GGV1 and GGV2).

### B MRST $`\lambda =0`$: spectral index

In our previous paper GGV2, we pointed out that an experimental measurement of the slope of the atmospheric lepton fluxes at energies where the prompt component dominates over the conventional one, might give information on the value of $`\lambda `$, the slope of the gluon PDF at small $`x`$. The best flux for this measurement is that of down-going muons, because the prompt neutrinos have first to convert into muons or electrons through a charged current interaction in order to be detectable in a neutrino telescope.

In this section and in the following two we consider the uncertainties in our method to determine $`\lambda `$. In this section we examine those coming from the charm production model, in Sect. V those related to the non linearity of our model, and in Sect. VI those coming from the unknown composition of the cosmic rays at high energies.

The slope of the fluxes or spectral index is $`\alpha _{\mathrm{}}(E_{\mathrm{}})=\mathrm{ln}\varphi _{\mathrm{}}(E_{\mathrm{}})/\mathrm{ln}E_{\mathrm{}}`$, with $`\mathrm{}=\mu ^\pm ,\nu _\mu +\overline{\nu _\mu }`$ or $`\nu _e+\overline{\nu _e}`$. In other words, the final lepton fluxes are

$$\varphi _{\mathrm{}}(E_{\mathrm{}})E_{\mathrm{}}^{\alpha _{\mathrm{}}(E_{\mathrm{}})}.$$
(21)
In GGV2 we found a simple linear dependence of $`\alpha _{\mathrm{}}`$ on $`\lambda `$, namely

$$\alpha _{\mathrm{}}(E_{\mathrm{}})=b_{\mathrm{}}(E_{\mathrm{}},\gamma ,\lambda )\lambda b_{\mathrm{}}(E_{\mathrm{}})\lambda ,$$
(22)
where $`b_{\mathrm{}}(E_{\mathrm{}})`$ is an energy dependent coefficient evaluated using our simulation with $`\lambda =0`$ and fixed $`\gamma `$. As argued in GGV2 (cfr. Eqs. (35) and (36) therein), the coefficient $`b_{\mathrm{}}(E_{\mathrm{}},\gamma ,\lambda )`$ depends mildly on $`\lambda `$ and can be well approximated by its value for $`\lambda =0`$ (see Sect. V). The coefficient $`b_{\mathrm{}}(E_{\mathrm{}},\gamma ,\lambda )`$ depends almost linearly on $`\gamma `$, the spectral index of the primary cosmic rays. We recall that the equivalent nucleon flux for primary cosmic rays is expressed as

$$\varphi _N(E)E^{\gamma 1}.$$
(23)
The linear dependence of $`b_{\mathrm{}}(E_{\mathrm{}},\gamma ,\lambda )`$ on $`\gamma `$ can be written as

$$b_{\mathrm{}}(E_{\mathrm{}},\gamma ,\lambda )=\overline{b}_{\mathrm{}}(E_{\mathrm{}},\gamma ,\lambda )+\gamma ,$$
(24)
where $`\overline{b}_{\mathrm{}}(E_{\mathrm{}},\gamma ,\lambda )`$ depends mildly on $`\lambda `$ and $`\gamma `$ We have included in $`\overline{b}_{\mathrm{}}`$ the $`+1`$ term coming from the $`1`$ in the exponent of Eq. (23)., as we will prove in Sect. V and Sect. VI, respectively.

Given $`b_{\mathrm{}}(E_{\mathrm{}})`$ from our model, an experimental measurement of $`\alpha _{\mathrm{}}`$ at energy $`E_{\mathrm{}}`$ would immediately give $`\lambda `$ corresponding to a value of $`x\mathrm{GeV}/E_{\mathrm{}}`$, as we discussed in GGV2. A measurement at $`E_{\mathrm{}}10^6`$ GeV = 1 PeV would give $`\lambda `$ at $`x10^6`$, a value of $`x`$ unattainable by present experiments.

For the time being we keep fixed the value of $`\gamma `$ ($`\gamma =1.7`$ below the knee, and $`\gamma =2.0`$ above the knee, as in GGV1 and GGV2); only in Sect. VI we will consider changes in the value of $`\gamma `$.

The feasibility of a measurement of $`\lambda `$ depends, therefore, on the uncertainties in $`b_{\mathrm{}}(E_{\mathrm{}})`$. Here we discuss those coming from the charm production model.

Fig. 4 shows the $`b_\mu `$ corresponding to the fluxes of Fig. 3 as functions of $`E_\mu `$. In the region of interest $`E_\mu 10^5`$$`10^6`$ GeV, the values of $`b_\mu `$ within each “band” decrease with increasing $`m_c`$ (and correspondingly smaller values of $`\mu _R`$), in Table I.

The spread of $`b_\mu `$ due to the choice of $`\mu _F`$ , $`\mu _R`$ and $`m_c`$ is $`\mathrm{\Delta }b_\mu 0.1`$, or $`\mathrm{\Delta }b_\mu /b_\mu 0.03`$, much smaller than the uncertainty $`\mathrm{\Delta }\varphi _\mu /\varphi _\mu 10`$ of the absolute flux in Fig. 3. If we would for some reason restrict ourselves to the $`\mu _F=2m_T`$ band, the uncertainty on $`b_\mu `$ would become even smaller, $`\mathrm{\Delta }b_\mu 0.03`$. We will refer to this error as $`\mathrm{\Delta }b_{par}`$ in the following, as it is related to the choice of parameters in the charm production model, and consider half of the spread in Fig. 4 to evaluate it. Therefore

$$\mathrm{\Delta }b_{par}0.05(0.015),$$
(25)
where the value in parenthesis corresponds to the $`\mu _F=2m_T`$ band.

### C MRST $`\lambda =\lambda `$(T)

So far we used $`\lambda =0`$ only. This case determines the uncertainty of the $`b_{\mathrm{}}(E_{\mathrm{}})`$ function which enters in the determination of $`\lambda `$ through the atmospheric muon fluxes.

Here we study an “intermediate” value of $`\lambda `$. We continue to use the MRST PDF, but with the value of $`\lambda =\lambda (T)`$ given by the slope of the lowest tabulated value of $`x`$ (see GGV2 for more explanations). This value depends on $`Q^2`$ and is about 0.2-0.3.

We repeat the same analysis of subsection IV A. However, for simplicity, we report the results for four selected sets of values for the parameters in Table I. The first set ($`m_c=1.1`$ GeV, $`\mu _F=0.5m_T`$, $`\mu _R=2.53m_T`$) gives a lower bound for the fluxes. The second set ($`m_c=1.4`$ GeV, $`\mu _F=2m_T`$, $`\mu _R=0.61m_T`$) gives an upper bound for the fluxes. The remaining two sets are cases in the $`\mu _F=2m_T`$ “band”.

The results are plotted in Fig. 5. The general features of Fig. 5 coincide with those of Fig. 3, except for an overall increase in all the fluxes due to the larger value of $`\lambda `$. The total spread of the fluxes given by the two limiting cases, as well as the spread within the narrower $`\mu _F=2m_T`$ band, are comparable to those found for $`\lambda =0`$. As in Fig. 3, our standard choice of parameters ($`m_c=1.25\mathrm{GeV},\mu _F=2m_T,\mu _R=1.0m_T`$) yields almost the highest flux.

We conclude that similar features would be obtained for other values of $`\lambda `$: our “standard choice” flux would be almost the highest in a band of uncertainty whose width is similar for all values of $`\lambda `$. The fluxes in the uncertainty band of Fig. 5 are consistent with older predictions (see GGV1 and references therein) and with the prediction by L. Pasquali et al. .

### D Other PDF’s

Another source of uncertainty for the final fluxes and spectral indices is the choice of the PDF set. As in GGV2, we consider here four recent sets of PDF’s: MRS R1-R2 , CTEQ 4M  and MRST , with the standard choice of parameters of Eqs. (13),(14),(15), (16),(17).

Figs. 6a and 6b show the muon fluxes (top panels) and spectral indices (bottom panels) for the two limiting cases of $`\lambda =0`$ (Fig. 6a) and $`\lambda =0.5`$ (Fig. 6b). In both cases the $`\mu `$ fluxes show at most a $`3050\%`$ variation depending on the PDF used. The uncertainty in the spectral indices for $`E_\mu 10^510^6`$ GeV is $`\mathrm{\Delta }b_\mu 0.02`$, or $`\mathrm{\Delta }b_\mu /b_\mu 0.01`$. This error will be denoted as $`\mathrm{\Delta }b_{PDF}`$ in the following, namely (again dividing the spread by $`2`$)

$$\mathrm{\Delta }b_{PDF}0.01.$$
(26)
These uncertainties, related to the PDF’s, are smaller that those due to the choices of mass scales (see Figs. 3-4). We conclude that, provided different PDF’s are calibrated in a similar way (i.e. same values of $`\mu _F`$, $`\mu _R`$ and $`m_c`$, chosen to fit the experimental data of our calibration), the final fluxes and spectral indices are very similar. The main source of uncertainty resides therefore in the choice of the mass parameters, rather than the adoption of a certain PDF set.

## V Determination of $`\lambda `$ with neutrino telescopes

In GGV2 we have given a detailed analysis of the dependence of the final fluxes and spectral indices on $`\lambda `$ for different PDF’s. In this section we show that the spread in our results due to $`\lambda `$ is larger than the one due to the choice of $`\mu _F`$, $`\mu _R`$, $`m_c`$ and of the PDF set, analyzed in the previous section. This is good news for the possibility of measuring $`\lambda `$, since the spread in $`\alpha _\mu `$, due to different $`\lambda `$’s, is the signal we want to detect, while the spread due to other factors constitutes the theoretical error of this measurement.

Figs. 7-10 show how the $`\mu `$ flux and its spectral index depend on $`\lambda `$. We used MRST with variable $`\lambda =0,0.1,0.2,0.3,0.4,0.5`$ and our standard choice of parameters ($`m_c=1.25\mathrm{GeV}`$, $`\mu _F=2m_T`$, $`\mu _R=1.0m_T`$).

Fig. 7 contains the differential muon fluxes. At the highest energies the $`\mu `$ fluxes are spread over almost two orders of magnitude. To each of the curves in this plot we need to assign a band of uncertainty of about one order of magnitude coming from the choice of the PDF and of the parameters $`m_c`$, $`\mu _F`$, and $`\mu _R`$ (see Fig. 3). Thus the curves become entirely superposed with each other. This makes it impossible to derive the value of $`\lambda `$ from an experimental measurement of the absolute level of the fluxes. However, the uncertainties in the spectral index of these prompt muons are much smaller and a determination of $`\lambda `$ becomes possible using the slope of the muon fluxes instead of their absolute level.

Fig. 8 shows the $`E^2`$-weighted integrated fluxes as functions of the muon energy. The slant lines indicate the number of particles traversing a km3 detector over a $`2\pi `$ sr solid angle. Even for the highest predicted fluxes, less than 1 particle per year will traverse the km3 detector for energies above $`10^8`$ GeV. Moreover, while prompt muons can be detected directly, prompt neutrinos have first to convert into charged leptons before being detected. The smallness of the charged current interaction effecting the conversion considerably lowers the detection rate of neutrinos. Therefore, the slope of the charm component of the atmospheric leptons can be studied in neutrino telescopes only using atmospheric muons coming from above the horizon, and only in a narrow range of energies, between a lower limit of $`E_\mu 10^510^6`$ GeV, above which the prompt component dominates over the conventional one, and an upper limit of $`E_\mu 10^710^8`$ GeV, above which the detection rates become negligible.

In practice, the spectral index of the prompt muon flux may be estimated by the difference of two integrated muon fluxes above two different energies, e.g. $`10^6`$ and $`10^7`$ GeV.

Figs. 9, 10 prove the validity in our model of Eq. (22), which is $`\alpha _{\mathrm{}}(E_{\mathrm{}})=b_{\mathrm{}}(E_{\mathrm{}})\lambda `$. In Fig. 9 we plot the spectral indices $`\alpha _{\mathrm{}}(E_{\mathrm{}})`$ for the different values of $`\lambda `$, both as directly calculated with our simulation (solid lines) and as $`b_{\mathrm{}}(E_{\mathrm{}})+\lambda `$ (dotted lines), where $`b_{\mathrm{}}(E_{\mathrm{}})`$ is $`\alpha _{\mathrm{}}`$ with $`\lambda =0`$. The two almost coincide, in the interval of interest, $`E_{\mathrm{}}10^6`$ GeV. Their difference, $`\alpha _{\mathrm{}}(E_{\mathrm{}})b_{\mathrm{}}(E_{\mathrm{}})+\lambda `$, given in Fig. 10, is small, about $`0.03`$ at $`E10^6`$ GeV. This difference stems from the mild dependence of $`b_{\mathrm{}}(E_{\mathrm{}})`$ on $`\lambda `$ and need to be added to the the other uncertainties evaluated in Sect. IV. We will refer to this error, due to the non linearity in $`\lambda `$ of Eq. (22), as

$$\mathrm{\Delta }b_{nonlin}0.015,$$
(27)
which again is evaluated dividing by $`2`$ the spread in Fig. 10.

We see in Fig. 9 that $`\mathrm{\Delta }\lambda \mathrm{\Delta }\alpha `$, therefore we would need an uncertainty in the spectral index of order $`0.1`$ to determine $`\lambda `$ with the same accuracy. We will show now that this is roughly the uncertainty related to our theoretical model.

In fact, we can finally estimate the total uncertainty in the determination of $`\lambda `$ coming from our theoretical model (that is, excluding the uncertainty due to the unknown composition of cosmic rays). We sum together the three spreads of $`b_{\mathrm{}}(E_{\mathrm{}})`$ previously calculated in Eqs. (25), (26) and (27), to obtain the final uncertainty We summed the errors linearly. Summing in quadrature would give $`(\mathrm{\Delta }\lambda )_{charm}(\mathrm{\Delta }b_\mu )_{charm}0.053(0.023)`$. from the charm production model,

$$(\mathrm{\Delta }\lambda )_{charm}(\mathrm{\Delta }b_\mu )_{charm}0.075(0.04),$$
(28)
where the number in parenthesis corresponds to fixing $`\mu _F=2m_T`$ in the charm model.

If the theoretical uncertainties so far presented would be the only ones affecting the determination of $`\lambda `$ through a measurement of the slope of the down-going muon flux, we could expect to get to know $`\lambda `$ with an uncertainty of about $`\mathrm{\Delta }\lambda 0.1`$. However, even excluding experimental uncertainties in the neutrino telescopes themselves, the uncertainty increases when our ignorance of the composition of the cosmic rays at high energy is taken into account, as we show in the following section.

## VI Uncertainty from cosmic ray composition

The final uncertainty we consider in the determination of $`\lambda `$ comes from the poorly known elemental composition of the high energy cosmic rays.

The spectral index of the cosmic rays $`\gamma `$ enters almost linearly in the slope of the atmospheric leptons. From Eqs. (22) and (24) we have

$$\alpha _{\mathrm{}}(E_{\mathrm{}})=\overline{b}_{\mathrm{}}(E_{\mathrm{}},\gamma ,\lambda )+\gamma \lambda .$$
(29)
So far we have kept $`\gamma `$ fixed, thus the uncertainty $`\mathrm{\Delta }b_\mu `$ calculated in Eq. (28) is actually an uncertainty in $`\overline{b}_\mu `$. We are going now to evaluate the uncertainty due to $`\gamma `$.

The non–linearity of Eq. (29) with respect to $`\gamma `$ is mild, as we have argued analytically in Sect. V of GGV2 and we show here using our numerical simulation. We have conducted a few trial runs of our simulation simply changing the values of $`\gamma `$ used for the primary flux. We recall from subsection IV B that in our model we used $`\gamma =1.7,2.0`$, respectively below and above the knee at $`E=510^6GeV`$. We have run our simulation changing these values of $`\gamma `$ by $`\pm 0.1,\pm 0.2`$ §§§Notice that these values of $`\gamma `$ are some of the most probable values (see Fig. 13)., both above and below the knee, to see the error produced when taking $`\overline{b}_{\mathrm{}}`$ computed at fixed $`\gamma `$ (our usual values) in Eq. (29) and thus leaving a pure linear dependence on $`\gamma `$. We used the MRST PDF, with $`\lambda =0`$, but similar results are obtained with other PDF’s and $`\lambda `$’s.

In Fig. 11a we plot the spectral index $`\alpha _{\mathrm{}}(E_{\mathrm{}})`$ for the different values of $`\gamma `$, both as directly calculated with our simulation (solid lines) and as $`\overline{b}_{\mathrm{}}(E_{\mathrm{}};\gamma =1.7,2.0;\lambda =0)\gamma `$ (dotted lines), i.e. using our standard values for $`\gamma `$ of the primary flux and adding an increment in $`\gamma `$ equal to $`\pm 0.1,\pm 0.2`$. In this way the “central value” of these spectral indices corresponds to the $`\lambda =0`$ case of Fig. 9. We can see that the dotted and the solid lines almost coincide, especially in the interval of interest for $`E_{\mathrm{}}10^510^6`$ GeV, proving the validity of Eq. (29). The uncertainty in $`\overline{b}_{\mathrm{}}`$ due to this non-linearity, that we call $`\mathrm{\Delta }\gamma _{nonlin}`$, evaluated in terms of the difference $`\alpha _{\mathrm{}}\overline{b}_{\mathrm{}}\gamma `$, is plotted in Fig. 11b and, in the energy range of interest, is

$$\mathrm{\Delta }\gamma _{nonlin}0.02.$$
(30)
We will now consider the error due to the poorly known elemental composition of the high energy cosmic rays. Concerning charm production, the relevant cosmic ray flux to be considered is the equivalent flux of nucleons impinging on the atmosphere. For a given cosmic ray flux, the equivalent flux of nucleons $`\varphi _{eq}(E_N)`$ depends in general on the composition of the cosmic rays. Nuclei with atomic number $`A`$ and energy $`E_A`$, coming with a flux $`\varphi _A(E_A)`$, contribute an amount $`A\varphi _A(AE_N)`$ to the equivalent flux of nucleons at energy $`E_N=E_A/A`$. So in total

$$\varphi _{eq}(E_N)=\underset{A}{}A\varphi _A(AE_N).$$
(31)
The uncertainties in the equivalent nucleon flux arise from the poorly known composition of cosmic rays in the energy range above the so-called knee, $`E_A10^6`$ GeV.

The actual $`\gamma `$ that enters into our proposed determination of $`\lambda `$ is the spectral index of the equivalent nucleon flux $`\gamma _{eq}`$, the equivalent cosmic rays spectral index for short. The equivalent nucleon flux is written as $`\varphi _{eq}E_N^{\gamma _{eq}1}`$, so that the spectral index $`\gamma _{eq}`$ is given by

$$\gamma _{eq}+1=\frac{E_N}{\varphi _{eq}}\frac{\varphi _{eq}}{E_N}=\frac{1}{\varphi _{eq}}\underset{A}{}A\varphi _A(\gamma _A+1),$$
(32)
where $`\gamma _A`$ is the spectral index of the component of atomic number $`A`$, i.e. $`\varphi _A(E_A)=k_AE_A^{\gamma _A1}`$.

We have calculated $`\varphi _{eq}`$ and $`\gamma _{eq}`$ using the experimental data of JACEE , CASA-MIA , HEGRA , and the data collected by Biermann et al, in Table 1 of Ref. , each with their respective compositions. Figs. 12 and 13 show the $`\varphi _{eq}`$ and the $`\gamma _{eq}`$ so obtained. Only the data of CASA-MIA  and HEGRA  reach energies $`E_N10^8`$ GeV, so we have not extended our analysis beyond $`10^8`$ GeV.

We have calculated the error band associated to $`\gamma _{eq}`$ in two different ways, because of the different parametrization of the composition used in Refs.  to . Refs.  give separate power law fits to the spectrum of each cosmic ray component,

$$\varphi _A(E_A)=k_AE_A^{\gamma _A1},$$
(33)
where the parameters $`k_A`$ and $`\gamma _A`$ have errors $`\mathrm{\Delta }k_A`$ and $`\mathrm{\Delta }\gamma _A`$. Standard propagation of errors gives, in this case,

$$\mathrm{\Delta }\varphi _{eq}=\left\{\underset{A}{}A^2\varphi _A^2\left[\left(\frac{\mathrm{\Delta }k_A}{k_A}\right)^2+\left(\mathrm{ln}(AE_N)\mathrm{\Delta }\gamma _A\right)^2\right]\right\}^{1/2}$$
(34)
and

$$\mathrm{\Delta }\gamma _{eq}=\left\{\underset{A}{}\frac{A^2\varphi _A^2}{\varphi _{eq}^2}\left[(\gamma _A\gamma _{eq})^2\left(\frac{\mathrm{\Delta }k_A}{k_A}\right)^2+\left[1(\gamma _A\gamma _{eq})\mathrm{ln}(AE_N)\right]^2\left(\mathrm{\Delta }\gamma _A\right)^2\right]\right\}^{1/2},$$
(35)
where $`\varphi _A`$ is evaluated at $`E_A=AE_N`$.

Refs. , give a power law fit to the total particle flux

$$\varphi (E_A)=kE_A^{\gamma 1}$$
(36)
and a composition ratio $`r_A(E_A)`$ in terms of which

$$\varphi _A(E_A)=r_A(E_A)\varphi (E_A).$$
(37)
These experiments distinguish only between a light and a heavy component. We assign atomic number 1 to the light component and 56 to the heavy one (which we call “iron”). Here $`k`$, $`\gamma `$, and $`r_A`$ have errors $`\mathrm{\Delta }k`$, $`\mathrm{\Delta }\gamma `$, and $`\mathrm{\Delta }r_A`$, respectively. The equivalent nucleon flux is still given by Eq. (31), while standard propagation of errors gives in this case

$$\mathrm{\Delta }\varphi _{eq}=\left\{\underset{A}{}A^2\varphi _A^2\left(\frac{\mathrm{\Delta }r_A}{r_A}\right)^2+\varphi _{eq}^2\left(\frac{\mathrm{\Delta }k}{k}\right)^2+\left[\underset{A}{}A\varphi _A\mathrm{ln}(AE_N)\mathrm{\Delta }\gamma _A\right]^2\right\}^{1/2},$$
(38)
We omit the much longer expression for $`\mathrm{\Delta }\gamma _{eq}`$. For simplicity, we have neglected the error coming from the energy dependence of $`r_A`$, which we expect to be much smaller than the others. In Fig. 12 we show the equivalent nucleon flux $`\varphi _{eq}`$. It is clear from the figure that the systematic uncertainties dominate, with spreads between different experiments of up to a factor of 4.

The uncertainties in the equivalent spectral index $`\gamma _{eq}`$ are smaller, as can be seen in Fig. 13, where only HEGRA and CASA-MIA extend to the energy region above the knee which is important to us.

We can consider, for example, an energy $`E_N10^7GeV`$, which is likely to determine the leptonic fluxes at around $`E_{\mathrm{}}10^6GeV`$, energy at which we would like to measure $`\lambda `$ through the spectral index (we recall from GGV2 that $`E_{\mathrm{}}0.1E_N`$).

At this energy $`E_N`$, from Fig. 13, we may take half the difference between the central values of the CASA-MIA and HEGRA data as an indication of the systematic uncertainty on $`\gamma _{eq}`$,

$$\mathrm{\Delta }\gamma _{syst}0.1.$$
(39)
Using the CASA-MIA data and the related error band, instead of the very spread HEGRA data, we can expect a reasonable statistical uncertainty

$$\mathrm{\Delta }\gamma _{stat}0.05.$$
(40)
Since $`\alpha _{\mathrm{}}`$ depends linearly on $`\gamma _{eq}`$ and $`\lambda `$, the same uncertainties apply to a determination of $`\lambda `$. The total uncertainty in the determination of $`\lambda `$ coming from the unknown composition of cosmic rays is now simply the sum of Eqs. (30), (39) and (40),

$$(\mathrm{\Delta }\lambda )_{comp}(\mathrm{\Delta }\gamma _{eq})_{comp}0.17,$$
(41)
if summing the errors linearly, or

$$(\mathrm{\Delta }\lambda )_{comp}(\mathrm{\Delta }\gamma _{eq})_{comp}0.11,$$
(42)
if we sum them in quadrature.

Finally, we can now combine all the uncertainties together, to compute the overall theoretical error in the determination of $`\lambda `$ with neutrino telescopes. From Eqs. (25), (26), (27), (30), (39), and (40) we obtain

$$\mathrm{\Delta }\lambda 0.25(0.21)$$
(43)
if summing errors linearly, or

$$\mathrm{\Delta }\lambda 0.13(0.11),$$
(44)
if summing in quadrature, where the numbers in parenthesis correspond to the $`\mu _F=2m_T`$ “band” in the charm model.

## VII Conclusions

We have examined in detail the possibility of determining the slope $`\lambda `$ of the gluon PDF, at momentum fraction $`x10^5`$, not reachable in laboratories, through the measurement in neutrino telescopes of the slope of down-going muon fluxes at $`E_\mu x^1GeV`$.

The method we are proposing may reasonably well reach $`x10^7`$, what would require $`10PeV`$ in muon energy. At this energy, there would still be $`50`$ events from charm if $`\lambda =0.5`$ and $`10`$ events if $`\lambda =0`$. But the best measurement could be done between $`100TeV`$ and $`1PeV`$ of muon energy, i.e. between $`x10^4`$ and $`x10^6`$. Present data do not go below $`x10^5`$ and the Large Hadron Collider (LHC) will not do any better. The reason is that the dominant values of $`x`$ in the production of a heavy particle of mass $`M`$ and rapidity $`y`$ are of order $`x[M\mathrm{exp}(\pm y)/\sqrt{s}]`$ (see for example ) where $`\sqrt{s}`$ is the center-of-mass energy of the hadron collision. Thus the smaller values of $`x`$ are obtained with the smaller $`M`$ and larger $`y`$ for fixed $`\sqrt{s}`$ ($`14TeV`$ at the LHC). Even if exhaustive studies of the possible minimum $`x`$ to be reached at the LHC have not yet been carried out , it is known that the experiments will explore the central rapidity region (the CMS and ATLAS detectors will cover $`y<0.9`$ only) and that bottom can be tagged, but most likely not charm Tagging is done by finding the point where the quark decays (called the vertex). The probability of decaying is exponential, and higher in the region close to the collision point. The only way to tell charm from bottom is by the distance from the collision to the vertex and bottom quark lives longer than charm. Thus, to detect charm with a good degree of confidence, one needs to select vertices close to the collision point. But in this region the vertices from bottom decay dominate, because the number of decay channels of the $`b`$ quark is five times larger than that of the $`c`$ quark. . This means that the lowest $`x`$ that LHC is expected to reach, assuming realistically that charm will not be tagged, is $`xm_b\mathrm{exp}(0.9)/\sqrt{s}=1.5\times 10^4`$. Therefore, the method proposed here may give information on the gluon PDF at $`x<10^5`$, a range not reachable in laboratory experiments in the near future.

To this end we studied the dependence of the leptonic fluxes and their slopes on $`\lambda `$. The slopes depend almost linearly on $`\lambda `$. We studied the uncertainties of the method we propose (excluding the experimental errors of the telescopes themselves). These come mainly from two sources: the free parameters of the NLO QCD calculation of charm production and the poorly known composition of cosmic rays at high energies.

We have seen that, for a fixed value of $`\lambda `$, the uncertainties give rise to an error band for the leptonic fluxes of almost one order of magnitude at the highest energies. This makes impossible a determination of $`\lambda `$ based solely on the absolute values of the fluxes, therefore we propose using the slopes of the fluxes. In particular we are proposing to use down-going muons, for energies $`E_\mu 100TeV`$, where prompt muons dominate over conventional ones, and not up-going neutrino-induced muons whose flux is orders of magnitude smaller. While an important contribution to up-going muons is expected from astrophysical neutrinos, there is no background for down-going atmospheric muons.

The overall theoretical error, from the charm production model, on the measurement of $`\lambda `$, is $`(\mathrm{\Delta }\lambda )_{charm}0.10`$. A comparable error, due to uncertainties in the cosmic ray composition, $`(\mathrm{\Delta }\lambda )_{comp}0.15`$, must be added, so that the overall error in the determination of $`\lambda `$ with neutrino telescopes is $`\mathrm{\Delta }\lambda 0.2`$.

These errors may be reduced by improving the experimental knowledge of the charm production cross sections and of the cosmic ray composition around and above the knee.

###### Acknowledgements.

The authors would like to thank M. Mangano, P. Nason, L. Rolandi and D. Treleani for helpful discussions. We also thank M. Mangano and P. Nason for the MNR program. This research was supported in part by the US Department of Energy under grant DE-FG03-91ER40662 Task C.

## FIGURE CAPTIONS

The ratio $`R=(\sigma _{NLO}\sigma _{LO})/(\sigma _{LO}\alpha _s\mathrm{ln}(s/m_c^2)/\pi )`$ is plotted as a function of the beam energy $`E`$, for the different values of $`\lambda `$ used with the MRST PDF.

Total cross sections for charm production $`\sigma _{c\overline{c}}`$, up to NLO, calculated with MRST ($`\lambda =0`$) and the values of $`m_c`$, $`\mu _F`$, $`\mu _R`$ of Table I, are compared with recent experimental values . For each “band” in the figures (i.e. for $`\mu _F=m_T/2`$, $`m_T`$ and $`2m_T`$) the cross sections increase with increasing $`m_c`$ (and correspondingly smaller values of $`\mu _R`$) in Table I (Fig. 2b is an enlargement of Fig. 2a).

Results for MRST $`\lambda =0`$. The $`E_\mu ^3`$-weighted vertical prompt fluxes, at NLO, are calculated using the values of $`m_c`$, $`\mu _F`$, $`\mu _R`$ of Table I and compared to the TIG  conventional and prompt fluxes. For each “band” in the figure (i.e. for $`\mu _F=m_T/2`$, $`m_T`$ and $`2m_T`$) the fluxes increase with increasing $`m_c`$ (and correspondingly smaller values of $`\mu _R`$) in Table I.

Spectral indices $`b_\mu `$ of the fluxes plotted in Fig. 3, for the MRST $`\lambda =0`$ case. For each “band” in the figure (i.e. for $`\mu _F=m_T/2`$, $`m_T`$ and $`2m_T`$) the spectral indices decrease with increasing $`m_c`$ (and correspondingly smaller values of $`\mu _R`$) in Table I.

Results for MRST $`\lambda =\lambda (T)`$. The $`E_\mu ^3`$-weighted vertical prompt fluxes, at NLO, are calculated using selected values of $`m_c`$, $`\mu _F`$, $`\mu _R`$ from Table I and compared to the TIG  conventional and prompt fluxes.

Results for MRS R1-R2, CTEQ 4M, MRST, for $`\lambda =0`$ (Fig. 6a) and $`\lambda =0.5`$ (Fig. 6b), with standard choice of parameters $`m_c`$, $`\mu _F`$, $`\mu _R`$. Top part: $`E_\mu ^3`$-weighted vertical prompt fluxes, at NLO. Bottom part: related spectral indices $`\alpha _\mu `$ (for the $`\lambda =0`$ case, $`\alpha _\mu =b_\mu `$).

Results for MRST $`\lambda =00.5`$ (solid lines). The $`E_\mu ^3`$-weighted vertical prompt fluxes, at NLO, are compared to the TIG  conventional and prompt fluxes (dashed lines).

Results for MRST $`\lambda =00.5`$ (solid lines). The $`E_\mu ^2`$-weighted integrated vertical prompt fluxes, at NLO, are compared to the number of particles traversing a km$`{}_{}{}^{3}2\pi `$ sr detector per year (dotted lines).

Results for MRST $`\lambda =00.5`$. The spectral indices $`\alpha _{\mathrm{}}(E_{\mathrm{}})`$ for the different values of $`\lambda `$, calculated directly by our simulation (solid lines) are compared to the corresponding terms $`b_{\mathrm{}}(E_{\mathrm{}})+\lambda `$ (dotted lines).

Results for MRST $`\lambda =00.5`$ (solid lines). The error of Eq. (22) is evaluated in terms of the difference $`\alpha _{\mathrm{}}(E_{\mathrm{}})b_{\mathrm{}}(E_{\mathrm{}})+\lambda `$.

Results for MRST $`\lambda =0`$ for different values of $`\gamma `$. a) The spectral indices $`\alpha _{\mathrm{}}(E_{\mathrm{}})`$ for the different values of $`\gamma `$, calculated directly by our simulation (solid lines) are compared to the corresponding terms $`\overline{b}_{\mathrm{}}(E_{\mathrm{}};\gamma =1.7,2.0;\lambda =0)\gamma `$, with increments in $`\gamma `$ equal to $`\pm 0.1,\pm 0.2`$ (dotted lines). The curves are labelled by the related value of $`\gamma `$ above the knee ($`\gamma =2.0`$ is our “standard value”). b) Uncertainty due to the non-linearity of Eq. (29), as the difference $`\alpha _{\mathrm{}}\overline{b}_{\mathrm{}}\gamma `$.

The $`E_N^3`$-weighted equivalent nucleon flux $`\varphi _{eq}(E_N)`$ is shown for different primary cosmic ray experiments . For each of these we plot the central value and the related error band.

The spectral index, $`\gamma _{eq}+1`$, for the equivalent nucleon fluxes of Fig. 12, is shown for different primary cosmic ray experiments . For each of these we plot the central value and the related error band.