text:
# Two Different Long-term Behaviors in Black-Hole Candidates: Evidence for Two Accretion Flows?

## 1 Introduction

Among the x-ray binaries which are identified as black-hole candidates by their mass functions, only three show persistent bright emission:  Cyg X–1, LMC X–3, and LMC X–1. At least two more non-transient Galactic objects are usually considered black-hole candidates by courtesy due to the similarity of their hard x-ray behavior to that of  Cyg X–1:  1E 1740.7–2942and  GRS 1758–258.  GX 339–4occupies a position intermediate between persistent and transient sources, being bright a significant fraction of the time but quiescent more often. The optical counterparts of the latter three systems have not been detected due to their distance, extinction by Galactic dust, and, in the case of  GX 339–4, dominance of the accretion-disk’s optical emission.

GRS 1758–258 almost certainly has a low or intermediate mass companion, and is therefore likely to be accreting by Roche-lobe overflow, in contrast to  Cyg X–1, which is thought to accrete the wind from its O supergiant companion. K-band observations of  GRS 1758–258 (Martí et al., 1998) show that the two brightest stars consistent with the Very Large Array (VLA) radio position appear to be a low-mass K-type red giant and a main-sequence F star; more massive stars, which might fuel wind accretion, are ruled out.

For  1E 1740.7–2942, Mirabel et al. (1991) suggested a companion later than B2, but their argument was indirect. The upper limit was on 6 cm radio flux; UV from an earlier companion would produce so much ionization in a positionally coincident molecular cloud that a compact HII region would have been seen with the VLA. This argument relies on  1E 1740.7–2942 being within the cloud and not behind or in front of it. HCO+ emission-line observations have been used to argue that it is indeed within the cloud (Philips & Lazio, 1995; Yan & Dalgarno, 1997), but the small amplitude of the 6.4 keV iron fluorescence line from the source has been used to argue that it is behind the cloud (Churazov, Gilfanov & Sunyaev, 1996) or in front of it (Sakano et al., 1999). The deepest K-band observations near  1E 1740.7–2942 were recently made with the Very Large Telescope (VLT) of the European Southern Observatory (Martí et al., 2000). Four stars of K magnitude 18.5-19.0 are about 1” from the VLA position, and each might represent a massive companion. Martí et al. (2000) found that the K-band 3$`\sigma `$ upper limit exactly at the VLA position (19.5) would constrain the companion to be a main sequence star later than about B8 V or a giant earlier than G5 III. The conclusion awaits improvement of the 1” uncertainty in the alignment of the VLA and VLT fields. Since  1E 1740.7–2942 shares the characteristic of bright, extended radio jets with  GRS 1758–258 and not  Cyg X–1, one might argue (independent of the evidence presented below) that it is more likely to be similar to  GRS 1758–258 overall, i.e. to be a Roche-lobe accretor.

Optical emission has been seen frequently from  GX 339–4, which shows less absorption than either  1E 1740.7–2942 or  GRS 1758–258, and is probably closer. The optical spectrum has always been consistent with an accretion disk, however, so there has been no spectral classification of the companion. Optical observations taken when the x-ray luminosity was very low (Callanan et al., 1992) imply an upper limit of (0.6 $`\pm `$ 0.2) L$`\mathrm{}`$ for the companion at the farthest estimated distance of 4 kpc (Cowley et al., 1987; Zdziarski et al., 1998) and one tenth of that at the closest estimated distance of 1.3 kpc (Predehl et al., 1991). Zdziarski et al. (1998) derived a somewhat higher extinction (Av = 3.4 - 4) than that taken by Callanan et al. (1992) from Ilovaisky et al. (1986) (Av = 2.2), which would increase the upper limit to about 3 L$`\mathrm{}`$. Even at this highest upper limit, the companion cannot be a high-mass star providing a significant wind.

The unique x-ray binary  Cyg X–3 has been interpreted as a black hole in a 4.8 hr orbit with a Wolf-Rayet star (Schmutz et al., 1996; van Kerkwijk et al., 1996), but this identification of both components has been questioned (Mitra, 1996, 1998). A Wolf-Rayet star could provide sufficient material in its wind to feed a bright x-ray binary system.

Cyg X–1,  1E 1740.7–2942,  GRS 1758–258, and  GX 339–4 spend most of their time in the spectral state called “hard” or “low”. This state, one of the two canonical states of black-hole candidates at moderate luminosity, is characterized by a hard power law (photon index 1.4 to 1.9) with an exponential rollover around 100 keV. The other canonical state, “soft” or “high”, is characterized by a softer power law with no exponential cutoff, plus a thermal component with a temperature on the order of 1 keV which can dominate the overall luminosity. There can be a gradual transition between these canonical states, now usually called an intermediate state (Méndez and van der Klis, 1997). It has recently been discovered, however, that at least one property of this state is not intermediate at all: Pottschmidt et al. (2000) found that the lags between rapid variations in soft and hard flux are far longer during state transitions of  Cyg X–1 than in either the soft or hard states. The “very high” state, characterized by high accretion rates and more complicated behaviors, is beyond the scope of this paper, never having been displayed by  Cyg X–1,  1E 1740.7–2942 or  GRS 1758–258.

One common picture of accretion in these systems is a standard thin accretion disk at large radii, which is truncated at some inner radius and replaced by a hot, more spherical flow (e.g., Shapiro et al., 1976). The thin disk produces thermal emission and the hot flow scatters it up into a power law. The transition from the low state to the high state is thought to occur when the accretion rate increases, with the result that the transition radius moves inward until the thin disk extends all the way down to the last stable orbit allowed by general relativity (e.g., Esin et al., 1998; Janiuk, Życki & Czerny, 2000).

In this paper we will demonstrate that this picture, although adequate to explain the behavior of  Cyg X–1, will have to be altered to address the more complicated evolution of luminosity and spectrum in other black-hole candidates.

## 2 Results

### 2.1 Observations

Cyg X–1,  Cyg X–3,  1E 1740.7–2942,  GRS 1758–258, and  GX 339–4 have been monitored by the Rossi X-Ray Timing Explorer (RXTE) since the start of science operations in 1996 February.  Cyg X–1,  Cyg X–3, and  GX 339–4 are bright enough to be monitored by the All-Sky Monitor (ASM) instrument, and we have used quick-look results provided by the ASM team.  1E 1740.7–2942 and  GRS 1758–258, which are fainter and in a crowded part of the sky, have been monitored by periodic visits with the highly-sensitive Proportional Counter Array (PCA) instrument (Jahoda et al., 1996).

The PCA monitoring observations, which last about 1500 s each, occurred monthly during 1996, weekly until early this year, and are now begin taken twice per week. There are data gaps from November to January every year due to a constraint against pointing the PCA near the Sun. In Main et al. (1999, hereafter Paper I) and Smith et al. (1997) we discussed the details of our observing strategy for  1E 1740.7–2942 and  GRS 1758–258, including offset-pointing to avoid nearby sources and subtraction of Galactic diffuse emission. The principal changes since Paper I was written are the new data acquired and a reanalysis of the whole data set using updated PCA response and background models (the new faint source model).

Figure 1 shows the count rate and best-fit power-law index (PLI) in photon space as a function of time for  1E 1740.7–2942 and  GRS 1758–258. The count rate is in the range 2.5 - 25 keV, uncorrected for interstellar absorption or instrument response. The PLI was fitted to the data in the same range.

Figure 2 shows ASM data for  Cyg X–1. To restrict our attention to the power-law component of the emission rather than the disk blackbody, we sum only ASM channels 2 (3.0-5.0 keV) and 3 (5.0-12.1 keV). To match the  1E 1740.7–2942 and  GRS 1758–258 data as closely as possible, we plot weekly accumulations of ASM data. Figure 2a shows the sum of the count rate in these two channels as a measure of the source flux. Figure 2b shows the ratio of these two channels.

We used weekly PCA pointings to  Cyg X–1 in 1998 to find the true PLI for comparison with contemporaneous ASM data. Figure 3 shows an expanded view of the ASM data with the PCA data superimposed. The count rates (Figure 3a) are for ASM channels 2+3 and for the PCA data restricted to the same range (3.0-12.1 keV). The PCA data show more variation not because of statistical errors, which are far smaller than the plotting symbols used, but because they represent a single (3000 s) snapshot instead of the average of a series of snapshots taken over a week, as is the case with the ASM data. This indicates a significant power at frequencies between about 1 and 300 $`\mu `$Hz, an interesting phenomenon beyond the scope of this paper.

Figure 3b compares the PCA PLI with the best-fit linear transform of the ASM channel ratio using this data set. This relation between the PLI and the ratio, $`R`$, of ASM channels 2 and 3 (Figure 2b) is

$$\mathrm{PLI}=1.499\mathrm{R}+0.698.$$
(1)
Figure 3c shows the PLI from the PCA data plotted against the ASM channel ratio for this data set, and the fit that gives equation 1. We do not expect that this relation applies for the softest data of Figure 2, where the disk blackbody emission could dominate channel 2, since this did not occur in 1998. The derived PLI, shown in Figure 2c, should therefore not be considered accurate when it is softer than about 2.1.

Figure 4 shows the ASM count rate and channel ratio for  GX 339–4. Since this source is much fainter than  Cyg X–1, and sometimes drops to undetectable levels, the data are poorer. The ASM coverage of any given point in the sky is irregular, so when the data are accumulated on a weekly basis as in Figure 2, there are periods where the statistics are so poor that a plot of the data is difficult to interpret. We therefore produced Figure 4 by accumulating sets of 50 consecutive ASM snapshots and averaging them together. One clear state transition can be see in 1998. Unfortunately, the transitions into and out of the soft state both occurred around January, when the RXTE pointing schedule, driven by a Sun-angle constraint, results annually in very poor coverage of the Galactic Center region by the ASM.

Finally, Figure 5 shows the ASM rate and ratio data for  Cyg X–3, with one-week accumulations as in Figure 2. The hard x-ray spectra of this source are more complicated than the others, including variable absorption and a very bright iron line, presumably due to the effects of a very dense wind from its companion (e.g., Nakamura et al., 1993). We have therefore not attempted to convert the channel 2 to channel 3 ratio into a power-law index, but we leave it as the measure of softness for this source.

### 2.2 Correlating Photon Flux and Spectral Index

There is a strong contrast in the behaviors shown in Figures 1 and 2. In  Cyg X–1, the curves of brightening and softening appear almost identical, while they have a very different time-dependence in  GRS 1758–258 and  1E 1740.7–2942. This phenomenon is the principal result of this paper. In order to make the relations between these quantities clearer in all three sources, we have derived Figures 6-9 from the raw data in Figures 1 and 2. These figures compare the PLI with both the photon flux and the time derivative of the photon flux for each source. In order to produce a set of standard comparisons, the count rates have had the instrumental responses and interstellar absorption removed and they are extrapolated to the range 2-100 keV using the fitted PLI. This range was chosen to represent the total photon flux in the power law part of the spectrum of each source. Above 100 keV the photon flux is small for both hard and soft state spectra, and below 100 keV the error introduced in the hard state by ignoring the exponential rolloff is also small. The flux derivatives were calculated at a given time by fitting a line to all the data within 15 days of that point and using the slope. For visual comparison, the PLI, flux, and flux derivative curves are all scaled to have a mean of zero and an rms of one.

The last panel of each figure shows the cross-correlation between the curves in each of the two panels above it. The peak correlations (positive or negative) and corresponding lags are shown in Table 1. A positive lag is defined as the PLI lagging behind the other quantity. To cross-correlate the unevenly sampled data, we interpolated between points on the PLI curve in order to get values at times corresponding to the shifted flux and flux derivative data.

The errors in the values of the correlation coefficient and the best-fit lag were estimated empirically. All points in the flux and PLI curves were randomly perturbed according to a Gaussian distribution with a width determined by the errors of the measurements. The analysis process, including calculation of the flux derivative and cross-correlation, was then repeated for the perturbed data set. This was done 1000 times with different random perturbations for each cross-correlation. The errors quoted on the correlation coefficients and the best-fit lags are the intervals that enclosed 95% of the results from the perturbed cases.

To make the results in Table 1 as mutually comparable as possible, we averaged the ASM data for  Cyg X–1 and  Cyg X–3 into time intervals centered on the times of each of the PCA pointings to  1E 1740.7–2942, so that the number of data points and irregularity of sampling is similar for all four sources. Figures 6, 7 and later figures for  Cyg X–1 and  Cyg X–3 use regular weekly binning of the data for display (as do Figures 2 and 5). Using these weekly bins or using daily bins, the cross-correlation results are still similar to Table 1. For  Cyg X–3 the ratio between ASM channels 2 and 3 is used as a surrogate for the PLI. For  Cyg X–1 the PLI and ASM channel ratio were shown to be linearly correlated (Figure 3), and therefore they yield the same curve when scaled to zero mean and unit rms.

Cyg X–1 and  Cyg X–3 (Figures 6 and 7) show an extremely good correlation of the PLI with the photon flux or count rate, but no clear correlation with its derivative.  1E 1740.7–2942 (Figure 8, Table 1) shows the opposite behavior: a very good correlation of the PLI with the derivative of the photon flux, and a poor correlation with the flux itself. The PLI leads the flux derivative significantly, i.e. the beginning of an episode of softening is a predictor of a drop in the photon flux.  GRS 1758–258 (Figure 9) shows behavior more like that of  1E 1740.7–2942 than that of  Cyg X–1, although the correlation with the flux derivative is not as good.

Qualitatively, the results of Figures 6-9 are the same whether the quantity compared to the PLI is the 2-100 keV photon flux as shown or the raw count rate in the instrument (PCA or ASM; see Table 1). When the energy flux from 2-100 keV is used in place of the photon flux,  1E 1740.7–2942 and  GRS 1758–258 still show a better correlation with the flux derivative than with the flux, but the correlations are weaker and the sign of the lag is reversed: the energy flux derivative leads the PLI (Figures 10 and 11). The corresponding plot for  Cyg X–1 using the energy flux looks similar to the photon flux plot (Figure 6). Because we don’t have a formula like equation 1 to convert the  Cyg X–3 ASM channel ratio to a PLI, we don’t calculate the energy flux for  Cyg X–3.

We don’t calculate a correlation for  GX 339–4 because there is only a single major episode in the plot shown in Figure 4: the rise to the soft state around January 1998 and the fall back to the hard state around December 1998. The first dotted line in Figure 4 shows the time at which the power law becomes soft. The second dotted line shows the time at which it hardens again. Note that at this time the flux has already dropped to a level considerably lower than that at which the power law first softened. Although this is a very small amount of data, qualitatively we might say that  GX 339–4 resembles  Cyg X–1 during the rise (softening correlated with flux) and  1E 1740.7–2942 and  GRS 1758–258 during the decline (a soft spectrum while the flux decreases). It has been noted that some black-hole x-ray nova outbursts have a hard-to-soft transition at high luminosity followed by a soft-to-hard transition at much lower luminosity (Miyamoto et al., 1995). In addition to citing the well-studied x-ray nova GS 1124-683 (Nova Muscae 1991), Miyamoto et al. found evidence for a non-simultaneous change of flux and PLI in  GX 339–4 using data taken by several instruments in 1988 and 1991. The total luminosity in the power law spectral component peaked as the PLI was making a transition from -2 to -3.

### 2.3 The most sudden changes

Although most of the large variations in count rate and PLI in  1E 1740.7–2942 and  GRS 1758–258 have been well sampled by our observing program (see Figure 1), each source has on at least one occasion changed too rapidly for us to follow the change in detail. Each of these events is worth special consideration, since the strongest tests of theoretical models for these sources may come from the most extreme variations.

Near the end of 2001 February, the power law flux from  GRS 1758–258 dropped by an order of magnitude between one pointing and the next, leaving behind a spectrum dominated by a weak blackbody component. Based on a model of two accretion flows (a sub-Keplerian halo and a thin disk), we predicted (Smith et al., 2001c) that mass input from the companion had ceased, and that the blackbody should decay smoothly over the following weeks. It did so, with a time constant of roughly 28 dy (Smith et al., 2001a). Both the sharp drop (due to the abrupt loss of the power-law component) and the subsequent decay of the soft flux are visible in the count rate plot of Figure 1. This decaying state was similar to the usual black hole “soft” or “high” state, both because it was dominated by thermal emission and because there was no detectable fast x-ray variability. A low-luminosity soft state for this source was previously inferred by Grebenev et al. (1997) from data taken between 1991 Fall and 1992 Spring: during this period, GRANAT found only upper limits above 3 keV, while ROSAT made a good spectral measurement from 1-2.4 keV. The combined data sets imply that the spectrum had to be both faint and soft. These results differ from the canonical “high/soft” state in that the soft spectrum appears at much lower luminosity than the hard state. The hard to soft transition is usually said to result from an increase in accretion rate and usually involves an increase in luminosity.

A single pointing on 2001 May 14 found that  1E 1740.7–2942 had suddenly jumped from one of softest PLIs ever seen from this source (2.23 $`\pm `$ 0.02) to one of the hardest (1.49 $`\pm `$ 0.02) (Smith et al., 2001b). As can be seen in Figure 12, the 2-100 keV (extrapolated) photon flux remained virtually unchanged while the spectrum hardened, so that the 2-100 keV energy flux briefly doubled, returning at the next pointing to its previous value. The raw PCA count rate (2.5-25 keV) increased slightly (this can be seen in Figure 1). If this one point was taken out of the data, the behavior in this period would be very similar to what we observe at other times: a softening preceeding a drop in the photon flux by about 10 dy (see above). We can be fairly certain that the May 14 result is not due to a second, very hard source turning on in the PCA field of view, because the hardening involved not just an increase in high-energy counts, but also a decrease in low-energy counts. Although this may be a rare sort of event, it might also be common but short-lived, and detected only once for that reason.

## 3 Discussion

### 3.1 Basic conclusions

Cyg X–1 has long been the prototype which theories of black-hole accretion have tried to model. The results above make it clear that there can be a richer variety of behavior even in other persistent black-hole binaries than is observed in  Cyg X–1. Specifically, any model in which a single variable (e.g. a single accretion rate) controls both the luminosity and the spectrum is clearly inadequate to explain the behavior of  1E 1740.7–2942 and  GRS 1758–258, even if it is extremely successful with  Cyg X–1.

The clearest contrast in behavior is between  Cyg X–1 and  1E 1740.7–2942: the PLI correlates with the photon flux in the former and its derivative in the latter.  Cyg X–3 is very much like  Cyg X–1 in this regard, and  GRS 1758–258 is more like  1E 1740.7–2942. More data are clearly needed to place  GX 339–4 firmly in either category. In $`\mathrm{\S }`$1, we reviewed the status of knowledge about the companion of each source. For  GRS 1758–258,  Cyg X–1, and  Cyg X–3, we can say that the sources with probable high-mass companions (wind accretors) have the PLI tracking with the flux, and the source with a probable low-mass companion has the PLI tracking the flux derivative. The situation is less clear for  GX 339–4 and  1E 1740.7–2942, in the former case because there are not enough data on the x-ray behavior, and in the latter case because there is neither a convincing measurement or upper limit for the companion.

If the correlation between companion type and x-ray behavior holds, the size of the accretion disk may be what distinguishes the two sorts of x-ray behavior. Wind accretors are expected to have a very small specific angular momentum in the accreting material, and therefore should form a very small disk. The complicated behavior in  1E 1740.7–2942 and  GRS 1758–258 might then be related to viscous propagation delays that are much smaller in  Cyg X–1. For a constant coefficient of kinematic viscosity, $`\nu `$, the viscous timescale in a standard thin disk $`t_{\mathrm{visc}}R^2/\nu `$ (Frank et al., 1992), so the accretion disk of  Cyg X–1 would only have to be an order of magnitude smaller than that of the LMXRBs to have a timescale $`<`$ 1 dy; i.e. it could still be over a thousand gravitational radii ($`R_g`$) in size. One recent model of wind-fed accretion (Beloborodov & Illarionov, 2001) produces a standard hard-state spectrum with a very quickly accreting, “inviscid” disk that is extremely small: only 14$`R_g`$. Compact objects which accrete via Roche-lobe overflow, however, have disks of a size on the order of the binary separation. These large disks have a viscous timescale (spiraling in of matter) which can be on the order of days to months (Frank et al., 1992), and could therefore be related to the time lags between correlated quantities in Figures 8 and 9. The $``$ 28-dy decay of  GRS 1758–258 in the soft state (Smith et al., 2001a) may be a direct measurement of the viscous timescale in that source.

Another significant difference between  Cyg X–1 on the one hand and  1E 1740.7–2942 and  GRS 1758–258 on the other is that the latter sources have bright radio lobes (Mirabel et al., 1992; Rodríguez, Mirabel, & Martí, 1992). This could be a consequence of disk size being related to jet collimation, of disruption of incipient jets by the companion wind in  Cyg X–1, or of some other phenomenon. If a substantial fraction of the hard x-ray emission in  1E 1740.7–2942 and  GRS 1758–258 occurred in the jets instead of near the disk, one might suggest that the time delays are due to the time required for information about changes in the accretion rate to propagate along the jet to the emitting region (e.g. a shock). The spectrum observed in  GRS 1758–258 as its emission shut off early in 2001 (Smith et al., 2001a) eliminates this possibility, however. That spectrum was dominated by a soft thermal component, brighter at first than any that was ever seen in the hard state. Since disk and jet spectral components would simply add, the thermal flux cannot be the jet emission since it is not usually seen in the hard state.

The hard-to-soft transition in  Cyg X–1 has sometimes been explained in the context of a central hot, advection-dominated flow with an outer thin disk, the transition radius determined by the accretion rate (e.g., Esin et al., 1998; Janiuk, Życki & Czerny, 2000). As the accretion rate increases, the central hot flow becomes unstable to collapse because it can cool itself more efficiently by bremsstrahlung. Here we will refer to this as a “static” soft state, since it can be reached by a gradual, quasi-static change in the accretion rate. In  1E 1740.7–2942 and  GRS 1758–258, the softest spectra observed are related to dropping photon fluxes (see Figures 8 & 9). Because the spectral state is related to the rate of change of luminosity, we refer to this as a “dynamical” soft state. Clearly it cannot be related to high accretion rates;  GRS 1758–258 showed a soft spectrum down to a few percent of its usual hard luminosity when it shut off (Smith et al., 2001a).

### 3.2 Relation to previous results

In Paper I, we had only the episodes of brightening and softening occurring in  1E 1740.7–2942 and  GRS 1758–258 in 1998 to look at. From these data alone, there was no way to tell whether the PLI was following the same profile as the photon flux but with a delay, or whether it instead followed the flux’s derivative (these options would of course be completely indistinguishable for sinusoidal variations). It is the newer data, including sharp drops in the count rate resulting in more significant softenings, which have decided the question.

We suggested two qualitative models in the Paper I for explaining the lack of simultaneity in the variations in flux and PLI. Both pictures involved the simultaneous presence of two accretion flows: a thin Keplerian disk and a hot, sub-Keplerian halo (e.g. an advection-dominated flow (Esin et al., 1998) or a shocked flow (Chakrabarti and Titarchuk, 1995)).

Our first picture followed a prediction of Chakrabarti and Titarchuk (1995): the sub-Keplerian flow might exist at all radii, not just within a transition radius, and if the mass accretion rate was boosted to both components at once, the halo would brighten at all radii almost instantaneously (near the free-fall timescale), while the inner regions of the thin disk would only brighten after the extra mass had wound its way in on the viscous timescale. We suggested that when this finally happened, the extra soft photons would cool the halo and the spectrum would soften.

Our second picture followed a suggestion by Mineshige (1996) related to the soft-to-hard state transition in the decaying phase of soft x-ray transients. He suggested that, for accretion rates where both thin-disk and advection-dominated flows are stable, the thin disk may tend to persist once it is established, i.e. the evaporation back to the advection-dominated flow is slow and inefficient. Thus, even if the brightening is quasi-static (i.e. the accretion rate changes slowly compared to the viscous timescale), there can still be a soft stage at the end of the outburst when the accretion rate is low again but the thin disk persists. Most disk analyses concentrate on finding stable solutions; the issue of transformations between stable states is difficult and much less well-explored.

Cyg X–1 provides a perfect test case to distinguish between these two pictures. Since it is expected to have a very short viscous timescale (see above), under the first picture (viscous delay), we would expect a much shorter delay. The second picture (the persistent disk) depends only on a hypothetical property of the inner disk, which should be roughly independent of the outer parts of the accretion flow. Thus, the simultaneity of brightening and softening in Figure 2 causes us to prefer the viscous delay picture (Chakrabarti and Titarchuk, 1995) over the slow-evaporation picture for the delay in  GRS 1758–258 and  1E 1740.7–2942.

Since we now have considerably more data, showing greater excursions in PLI, a deeper examination of the two-flow picture (Chakrabarti and Titarchuk, 1995) is warranted.

### 3.3 Correlating photon flux and derived halo depth

The hard component in black-hole-candidate spectra has long been attributed to inverse Compton scattering of soft disk photons by energetic electrons, usually but not always thermal. These hot electrons have been hypothesized to reside in several places, often a hot, central region of the disk or a halo or corona above the disk. The PLI of the Comptonized component is a function of the temperature and optical depth of the Comptonizing plasma. Pietrini & Krolik (1995) have modeled coronae above the disk and found that, for accretion rates much less than Eddington, the halo temperature is nearly constant at a value of about kT = 100 keV even as the accretion rate varies by up to four orders of magnitude. A similar result was found by Chakrabarti and Titarchuk (1995), but the hot region was a shocked plasma, originating from a sub-Keplerian flow, and occupying the radially innermost part of the accretion flow. The constant temperature holds as long as the dominant cooling mechanism for the halo is the Comptonization of soft photons and the accretion is significantly sub-Eddington. The highest luminosities derived from the  1E 1740.7–2942 and  GRS 1758–258 data of Figure 1 are 9.8 and 5.4 $`\times 10^9`$ erg cm-2 s-1 (extrapolated to 2-100 keV and with interstellar absorption removed). For a Galactic Center distance of 8.5 kpc, this gives 8.5 and 4.7 $`\times 10^{37}`$ erg s-1, or about 7% and 4% of the Eddington luminosity for a 10 M black hole, so these objects could plausibly be in the constant-temperature regime, in which the PLI depends only on the optical depth $`\tau `$ of the scattering region.

In unsaturated Comptonization, the PLI (Shapiro et al., 1976) is

$$\gamma =\left(\frac{9}{4}+\frac{4}{y(1+f)}\right)^{1/2}\frac{1}{2}$$
(2)
where the spectrum in photons cm-2 s-1 keV$`{}_{}{}^{1}E^\gamma `$. The Comptonization parameter

$$y=\frac{4kT_e}{m_ec^2}\mathrm{Max}(\tau ,\tau ^2)$$
(3)
and $`f`$ is a function only of the halo electron temperature $`T_e`$:

$$f=\mathrm{\hspace{0.33em}2.5}\theta +1.875\theta ^2(1\theta ),\theta kT_e/m_ec^2$$
(4)
which we are taking as constant. Then, at $`kT_e=`$100 keV,

$$\tau =\mathrm{\hspace{0.33em}3.3}\left[(\gamma +\frac{1}{2})^2\frac{9}{4}\right]^1\mathrm{or}\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}3.3}^{1/2}\left[(\gamma +\frac{1}{2})^2\frac{9}{4}\right]^{1/2},$$
(5)
whichever is less. With this assumption, $`\tau `$ for  1E 1740.7–2942 and  GRS 1758–258 remains close to 1, varying only from 0.44 to 1.39. Since $`\tau `$ increases monotonically as $`\gamma `$ decreases, its correlations with photon flux, energy flux, and their derivatives are of the opposite sign and similar magnitude (Figures 13-17 and Table 1).

For an accreting halo flow with an unchanging geometry and constant infall speed, $`\tau `$ is proportional to the accretion rate in this flow, $`\dot{M}_\mathrm{H}`$. Thus, if the geometry, infall speed, and temperature of the halo are all constant, the PLI (via $`\tau `$) is a surrogate for $`\dot{M}_\mathrm{H}`$. If we could additionally postulate that the photon flux $`F`$ is a surrogate for the flow in the thin disk, $`\dot{M}_\mathrm{D}`$, and that the hot halo is an independent flow, then Figures 13-15 could be interpreted in simple and compelling ways. For  1E 1740.7–2942 and  GRS 1758–258, where we have a good correlation in the second panel, we would have

$$\tau (t+\delta )\overline{\tau }\frac{d}{dt}F(t)\dot{M}_\mathrm{H}(t+\delta )\overline{\dot{M}_\mathrm{H}}\frac{d}{dt}\dot{M}_\mathrm{D}(t),$$
(6)
where $`\delta <0`$ is the lag shown in Table 1. The mean value is subtracted from $`\tau `$ because of the offset to zero average imposed on all the curves in these figures. This offset has little effect on $`dF/dt`$, since it already averages nearly to zero, there being no major long-term trends in the count rates.

Equation 6 suggests that the flow in the thin disk at any time is an integration of the halo flow in the recent past. This would be expected if there were two independent accretion flows: a thin disk (with its slow time constant) and the halo (nearly in free fall), and if they were supplied with proportional amounts of matter at large radii at any given time. The possibility of complicated variations of PLI and luminosity as a result of independent disk and halo flows was predicted by Chakrabarti and Titarchuk (1995). Integration and delay are both natural consequences of viscous propagation in a thin disk (Frank et al., 1992).

Under these assumptions (constant halo temperature and $`F`$ proportional to disk accretion rate), the observed behavior of  Cyg X–1, with a strong inverse correlation in the top panel of Figure 13, could result from

$$\dot{M}_\mathrm{H}(t)+\dot{M}_\mathrm{D}(t)\mathrm{constant}.$$
(7)
In other words, the state transitions could be due to the mass input switching back and forth between going mostly into the halo flow and mostly into the disk. This is an appealing picture of  Cyg X–1 for two reasons. First, the bolometric luminosity is, in fact, very nearly constant during state transitions of  Cyg X–1 (Zhang et al., 1997), which is not the case for other black-hole candidates. Second, models which attribute the state changes to a change in a single accretion rate (e.g., Esin et al., 1998; Janiuk, Życki & Czerny, 2000) require the average accretion rate to be in the narrow range of values near where the hard-to-soft transition naturally occurs. If the state transitions are instead due to the swapping of the accretion flow between two modes, they could occur over a wide range of total accretion rates.

The principal problem with these simple and compelling pictures is the use of $`F`$ as a surrogate for the thin disk flow. If the thin disk were a point source of soft photons at the center of a large, optically thick Comptonizing cloud, then every photon observed in the power law would have its origin in the thin disk, and (nearly) every photon emitted by the thin disk would end up in the power law. Such a halo has been postulated by Hua, Kazanas & Cui (1999) to explain the frequency dependence of lags between soft and hard photons at sub-second timescales in black-hole candidates. If the thin disk were a blackbody or a set of concentric blackbodies at different temperatures, then $`F`$ (first from the thin disk and, after Comptonization, in the power law) would go as $`\dot{M}_\mathrm{D}^{4/3}`$, not very different from $`\dot{M}_\mathrm{D}`$. For a halo which is close to the disk and tenuous, however, as is sometimes the case in recent models, $`\tau <<1`$, most disk photons are not Comptonized, and $`F`$ in the power law can be at least as strong a function of $`\dot{M}_\mathrm{H}`$ as $`\dot{M}_\mathrm{D}`$. A thin halo near the disk also re-heats the disk with some of the Comptonized flux, further complicating the issue. That influence would be lessened if the halo flow were geometrically large, so that most of the Comptonized flux does not intersect the thin disk, or if it had relativistic bulk motion away from the disk, so that the Comptonized photons moved preferentially outward (Beloborodov, 1999).

Finally, viscous action in the disk should not only delay but also smooth out variations in the input to the accretion disk. If the $`F`$ were a good surrogate for the disk accretion flow, it would therefore have less overall variability than the halo flow (represented for the sake of this discussion by $`\tau `$). This is not the case, however. For $`\tau 1`$, the first exponent in equation 5 is applicable. In this case, the rms variability in $`\tau `$ and in $`F`$ are remarkably similar. The rms variabilities in $`F`$ and $`\tau `$ are 19% and 18% respectively in  1E 1740.7–2942 and 22% and 21% in  GRS 1758–258. This can’t be seen in Figures 14-15 (top panels) because of the scaling used, which automatically sets the rms to 1, but the result is that for $`\tau 1`$ the upper panels in these figures would look the same even without the re-scaling. For $`\tau >1`$, the second exponent in equation 5 applies and the problem is aggravated: the rms variabilites in $`\tau `$ are only 9.4% in  1E 1740.7–2942 and 11% in  GRS 1758–258, where we would expect them to be higher than the variabilities in $`F`$.

On the other hand, the event of 2001 May 14 in  1E 1740.7–2942 (Figure 12) involves both the PLI and energy flux undergoing a large change while the photon flux remains unchanged. This shows that $`F`$ can sometimes have the slowest response to sudden changes, an expected characteristic of the inner disk accretion rate.

Our purpose, however, is not to argue that a very large, optically thick halo with a uniform temperature actually exists. Rather, it is to stimulate more realistic theoretical work, and to point out that any model of accretion in black-hole binaries must have at least the potential to explain the behavior epitomized by Figure 8. We have used a model of two independent flows and the postulated equivalence in equation 6 as a starting point. In a realistic two-flow model, the observables $`F`$ and $`\gamma `$ would both be functions of both $`\dot{M}_\mathrm{D}`$ and $`\dot{M}_\mathrm{H}`$. If the two-flow idea is valid, making these relations more physical would result in maintaining the time lag between $`\dot{M}_\mathrm{H}`$ and $`\dot{M}_\mathrm{D}`$ while making the rms variability of the former greater than the latter.

The disk and halo flows need not have an independent origin as they do in Chakrabarti and Titarchuk (1995). For example, Watanabe & Fukue (1996) discussed how an initially Keplerian corona, produced by evaporation from the disk, could lose angular momentum and accrete rapidly due to radiation drag from a central luminous source (i.e. a neutron star). Greater luminosity from the outer Keplerian thin disk would stabilize against this effect. For a black-hole system, we take the central radiation source as being the innermost region of the disk itself. If the mass input to the thin disk dropped suddenly, the inner disk radiation would at first dominate over the outer disk radiation since the drop would not yet have propagated inward. Thus the corona would accrete more quickly than normal. If the corona is produced at a constant rate by evaporation, the result would be a thinner corona and a softer power-law. Once the drop in accretion rate had propagated to the innermost regions, the former disk/corona ratio would re-establish itself and the spectrum would harden again at a lower luminosity. This reproduces the qualitative behavior in Figures 8 and 9: a softening precedes a drop in flux. What this picture has in common with Chakrabarti and Titarchuk (1995) is the viscous delay in the thin disk as the cause of the complication in the observed behavior.

## 4 Summary

We have found that black-hole candidates  GRS 1758–258 and  1E 1740.7–2942 display long-term behavior very different from the well-known variation of  Cyg X–1. The primary result, that the PLI varies as the negative derivative of the power-law flux instead of directly with it, is at least suggestive of a model with both slow (disk) and fast (halo) accretion flows.

Our observations come at a time when there is mounting evidence that even the outer parts of accretion disks contain a vertically extended flow as well as a standard thin disk. Cottam et al. (2001) observed the low-mass x-ray binary EXO 0748–67 with the Reflection Grating Spectrometer on the XMM-Newton observatory. Although this is an eclipsing binary, they found very bright x-ray emission lines which do not vary with the variation of the continuum from the central source. From this lack of variability and the details of the absorption and emission in the spectrum, they concluded that there must be a significant flow at large radii with a height much greater than that predicted by hydrostatic equilibrium of a thin disk. Smale et al. (2000) observed another dipping binary, X1624–490, with RXTE. They found that the Comptonization region which produces the hard x-ray continuum is never completely occulted during dips, while the blackbody from the central source is. From the broad profile of the dips, they deduced that the hard-x-ray-emitting region was $`5\times 10^{10}`$ cm in radius or larger, within a factor of two of the full presumed size of the thin accretion disk based on the 21-hour orbit of the system.

Kaaret et al. (1998) described certain correlations in the atoll sources (neutron-star binaries) 4U 0614+091 and 4U 1608–52 that might also be due to independent disk and halo flows. In those sources, the frequency of the kilohertz quasi-periodic oscillations (QPOs) correlates well with the power-law index (Kaaret et al., 1998) and, for 4U 0614+091, with the disk blackbody flux (Ford et al., 1997). Neither these sources, nor the kilohertz QPO sources in general, show a long-term correlation between QPO frequency and luminosity. Kaaret et al. (1998) suggested that the thin-disk accretion rate may control the blackbody luminosity, the QPO frequency, and (via Compton cooling of the halo) the power-law index, while the overall luminosity (dominated by the power-law) is controlled by the halo accretion rate.

Our secondary conclusion is that there can be two causes of spectral softening in black-hole binaries: the well-modeled static soft state, in which the accretion rate becomes high enough to force the innermost part of the halo into the thin disk configuration, and a newly-discovered dynamical soft state, occurring whenever the photon flux is dropping. In the context of a picture of two independent accretion flows, the dynamical soft state would occur when the inner thin disk has yet to respond to a drop in accretion which has already depleted the halo.

Because persistent black hole candidates are rare, there are only a few more immediately available tests of the ideas developed here. LMC X–3 is one source we have not yet examined. It is too faint for the RXTE ASM data to be useful spectrally, and the PCA monitoring observations so far, although they have revealed that the source does undergo soft/hard transitions (Wilms et al., 2001), have been too few to provide evidence for or against viscous delay. With a main sequence companion of only modestly high mass (B3 V), it is probably a disk accretor, and may therefore be a good candidate for delays.  GX 339–4, which we expect to show delays, may provide more conclusive data than the hint visible in Figure 4. A new PCA monitoring campaign for this source began in 2001 March. Finally, the XMM-Newton x-ray telescope has sufficient sensitivity to measure spectra of sources as faint as  1E 1740.7–2942 and  GRS 1758–258 in M31 with a reasonable exposure time ($`<`$ 1 dy). A monitoring campaign of the M31 bulge region could reveal a number of new sources with behavior similar to  1E 1740.7–2942 and  GRS 1758–258. If this behavior is indeed due to viscous delay in the disk, this campaign would produce a census of high-mass versus low-mass binaries using viscous delay as a distinguishing diagnostic.

In future work we will consider the connection between the behavior exhibited by  1E 1740.7–2942 and  GRS 1758–258 and that observed in the x-ray novae or soft x-ray transients. The transient systems, exemplified by Nova Muscae 1991 (GS 1124–683), tend to brighten in their hard state, change to the soft state near their peak, and decay in the soft state, returning to the hard state only at much lower luminosities (Miyamoto et al., 1995).  GX 339–4, much more variable than  1E 1740.7–2942 and  GRS 1758–258 but more persistent than the transients, may serve to clarify the connection.

The authors thank two anonymous referees for many useful suggestions and corrections. This work was funded in part by NASA grant NAG5-7265.


format: text

text:
# Atmospheric neutrino flux supported by recent muon experiments

## 1 Introduction

The atmospheric neutrino (AN) flux is one of the main tools for studying neutrino oscillations, neutrino decay and nonstandard neutrino-matter interactions. Furthermore it represents an unavoidable background for key experiments in astroparticle physics, as those searching for proton decay, $`n\overline{n}`$ transitions in nuclei, etc. Therefore, the accuracy of theoretical predictions of the AN flux ultimately affects the interpretation of many experimental data.

After several studies performed in the frameworks of one-dimensional (1D)  and three-dimensional (3D)  cascade models, the physics of neutrino production in the sub- and multi-GeV energy ranges is fairly well understood. However, the validity of the AN flux calculations is still controversial. This is mainly related to the uncertainties in the required input data (inclusive cross sections for particle production, primary cosmic ray spectrum and composition, etc.). Moreover, the problem is too complex to be solved without many simplifications (explicit or implicit) and these also are essential sources of uncertainties in the predicted AN fluxes . A popular opinion is that the accuracy of the present-day calculations of the low-energy AN flux is about 20–30%. Sometimes this belief of obscure origin is treated as the uncertainty in the normalization factor. But, as a matter of fact, the uncertainties in the input data are so large that the actual accuracy of the AN flux turns out to be much worse. It is also essential that the uncertainties are very energy dependent and therefore they cannot be reduced to a simple renormalization.

In this work we present a new 1D calculation of low and intermediate energy muon and neutrino fluxes, based on up-to-date data on primary cosmic-ray flux and hadronic interactions. Our calculation is based on a kinetic approach. In order to check the validity of our description of hadronic interactions and shower development, we perform a comparison of the predicted atmospheric muon fluxes with the data from several recent balloon-borne experiments.

We believe that an accurate 1D AN flux calculation is still helpful, despite the results of 3D calculations are now available. Indeed, the main lessons we got from the recent 3D analyses  may be summarized as follows.

* 3D effects lead to a strong modification of zenith and azimuth angle distributions of neutrinos at energies $`E_\nu 0.5`$ GeV. On the other hand, the corresponding effect for the neutrino-induced events in underground detectors is expected to be modest and scarcely measurable, because their distributions are close to isotropic out of a very big mean scattering angle in quasielastic neutrino-nucleus collisions at low energies.
* Even at very low energies the 3D corrections to the $`4\pi `$ and even $`2\pi `$ averaged AN fluxes are relatively small. At least, they are comparable with the uncertainties in the input data, including the uncertainties in the cross sections for exclusive neutrino-nucleus interactions needed for simulating the neutrino-induced events in the underground detectors.
* Above $`12`$ GeV the 3D effects become negligible for the neutrino fluxes averaged over the azimuth angles at any zenith angle.

On that ground, we expect the calculation offered for consideration to be valid without restrictions at $`E_\nu 1`$ GeV and we consider it as a basis for future 3D generalization of our approach to the AN problem.

Furthermore, our 1D code is useful for testing 1D and 3D codes with a method completely different from Monte Carlo. It can be used to perform a systematic analysis of uncertainties in the AN flux calculations and to adjust some poorly known input parameters, using the data on cosmic-ray secondaries. In this work, in order to evaluate quantitatively the uncertainties related to the particle interactions, we perform a comparison of the results obtained with different models for meson production in nucleon-nucleus and nucleus-nucleus collisions.

## 2 The CORT code

The present work is based on an updated code CORT (Cosmic-Origin Radiation Transport), previously developed in  and used in  (see also ). Like the earlier version, the current Fortran 90 code implements a numerical solution of a system of 1D kinetic equations describing the propagation of nuclei, nucleons, light mesons, muons and (anti)neutrinos of low and intermediate energies through a spherical, nonisothermal atmosphere. It takes into account solar modulation and geomagnetic effects, energy loss of charged particles, muon polarization and depolarization effect. The exact kinematics is used in description of particle interactions and decays. Let us sketch some most distinctive aspects.

In order to evaluate geomagnetic effects and to take into account the anisotropy of the primary cosmic-ray flux in the vicinity of the Earth, we use the method of ref.  and detailed maps of the effective vertical cutoff rigidities from ref. .111The penumbral effects are included in the definition of the effective cutoff. The maps are corrected for the geomagnetic pole drift and compared with the later results reviewed by Smart and Shea  and with the recent data on the proton flux in near earth orbit obtained with the Alpha Magnetic Spectrometer (AMS) . The interpolation between the reference points of the maps is performed by means of two-dimensional local B-spline. The Quenby-Wenk relation (see e.g. ), re-normalized to the vertical cutoffs, is applied for evaluating the effective cutoffs for oblique directions. More sophisticated effects, like the short-period variations of the geomagnetic field, Forbush decrease, re-entrant cosmic-ray albedo contribution, etc., are neglected. We also neglect the geomagnetic bending of the trajectories of charged secondaries and multiple scattering effects. Validity of our treatment of secondary nucleons and nuclei was confirmed using all available data on proton spectra in the atmosphere (see ).

The meteorological effects are included using the Dorman model of the atmosphere  which assumes an isothermal stratosphere and constant gradient of temperature (as a function of depths) below the tropopause. Ionization, radiative and photonuclear muon energy losses are treated as continuous processes. This approximation is quite tolerable for atmospheric depths $`h2\times 10^3`$ g/cm2 at all energies of interest . Propagation of $`\mu ^+`$ and $`\mu ^{}`$ originating from every source (pion or kaon decay) is described by separate kinetic equations for muons with definite polarization at production. These equations automatically account for muon depolarization through the energy loss (but not through the Coulomb scattering).

In the new version, essentially all approximate formulas used in  are replaced (or duplicated for a control) with the more accurate numerical equivalents.

## 3 Primary cosmic ray spectrum and composition

In the present calculations, the nuclear component of primary cosmic rays is broken up into 5 principal groups: H, He, CNO, Ne-S and Fe with average atomic masses $`A`$ of 1, 4, 15, 27 and 56, respectively. We do not take into account the isotopic composition of the primary nuclei and assume $`Z=A/2`$ for $`A>1`$, since the expected effect on the secondary lepton fluxes is estimated to be small with respect to present-day experimental uncertainties in the absolute cosmic-ray flux and chemical composition.222However, a more rigorous analysis must allow at least for small fractions of primary 2H and 3He affecting the $`\mu ^+/\mu ^{}`$ and $`\nu /\overline{\nu }`$ ratios.

We use the following 9-parameter model for the differential energy spectra of every nuclear group at the top of the atmosphere

$$\frac{F_\mathrm{A}(E)}{E}=\{\begin{array}{ccc}F_\mathrm{A}^1ϵ_\mathrm{A}^{\mathrm{\Gamma }_\mathrm{A}(ϵ_\mathrm{A})}& \text{at}& E<W_\mathrm{A},\\ F_\mathrm{A}^2ϵ_\mathrm{A}^{\gamma _\mathrm{A}}& \text{at}& E>W_\mathrm{A},\end{array}$$
(1)
with

$$\mathrm{\Gamma }_\mathrm{A}(ϵ_\mathrm{A})=\underset{k=0}{\overset{3}{}}\gamma _\mathrm{A}^{(k)}\mathrm{ln}^kϵ_\mathrm{A},ϵ_\mathrm{A}=r_\mathrm{A}E/W_\mathrm{A}.$$
(2)
Here $`E`$ is the total energy in GeV/nucleon. The condition for continuity of the spectrum and its first derivative in the point $`E=W_\mathrm{A}`$ yields

$$\mathrm{ln}\left(F_\mathrm{A}^2/F_\mathrm{A}^1\right)=\underset{k=1}{\overset{3}{}}k\gamma _\mathrm{A}^{(k)}\mathrm{ln}^{k+1}r_\mathrm{A},\gamma _\mathrm{A}=\underset{k=0}{\overset{3}{}}(k+1)\gamma _\mathrm{A}^{(k)}\mathrm{ln}^kr_\mathrm{A}.$$

Thus, we remain with 7 independent parameters for each group which we determine by fitting the experimental data.

For the hydrogen and helium groups at $`E<120`$ GeV/nucleon we chose the data of the balloon-borne experiment BESS obtained by a flight in 1998  (a period close to a minimum of solar activity). For higher energies (but below the knee) we use data by a series of twelve balloon flights of JACEE  and the result of an analysis by Wiebel-Sooth et al.  based upon a very representative compilation of world data on primaries and solid theoretical considerations. Since the fits obtained in  turned out to be very close to the combined result of the JACEE experiments, we call our model “BESS+JACEE” fit. A maximum likelihood analysis with $`W_1=W_4=2`$ TeV/nucleon gives the following values for the remaining parameters in eqs. (1) and (2):

$`F_1^1`$ $`=5.7216\times 10^5,`$ $`F_4^1`$ $`=1.6500\times 10^7,`$
$`F_1^2`$ $`=5.8511\times 10^5,`$ $`F_4^2`$ $`=1.7239\times 10^7,`$
$`r_1`$ $`=1.5952,`$ $`r_4`$ $`=0.4783,`$
$`\gamma _1`$ $`=2.800,`$ $`\gamma _4`$ $`=2.655,`$
$`\gamma _1^{(0)}`$ $`=2.7119,`$ $`\gamma _4^{(0)}`$ $`=2.7881,`$
$`\gamma _1^{(1)}`$ $`=6.9771\times 10^2,`$ $`\gamma _4^{(1)}`$ $`=1.1953\times 10^1,`$
$`\gamma _1^{(2)}`$ $`=3.2396\times 10^2,`$ $`\gamma _4^{(2)}`$ $`=2.8694\times 10^2,`$
$`\gamma _1^{(3)}`$ $`=3.8561\times 10^3,`$ $`\gamma _4^{(3)}`$ $`=2.0829\times 10^3,`$
where $`F_\mathrm{A}^i`$ are in m-2s-1sr-1(GeV/nucleon)-1. Clearly the many digits in the above parameters are by no means indicating the accuracy of the data, but they are necessary to avoid loss of accuracy in numerical calculations.

Figure 1 shows a comparison between the BESS+JACEE fit, the data from , the fit from  (shaded areas) and several other experiments , performed in different periods of solar activity. The filled/shaded areas in fig. 1 represent the power-type parametrizations derived in from the original data. The legend indicates the publication dates and (when relevant) the years of measurements.

One sees that the BESS+JACEE fit is in excellent agreement with the recent AMS data on the proton flux . On the other hand it slightly exceeds the AMS helium flux  (with average discrepancy of about 15%). Our choice of the BESS 98 data was in particular conditioned by the fact that the power-law extrapolation of the high-energy tail of the AMS helium rigidity spectrum ($`R^{(2.74\pm 0.02)}`$)  distinctly worse joins the world-averaged fit ($`E^{(2.64\pm 0.02)}`$)  and even the JACEE 1-12 fit ($`E^{\left(2.68_{0.06}^{+0.04}\right)}`$) .333This is all the more true for the remaining data presented in fig. 1.

We assume that the spectra of the remaining three nuclear groups are similar to the helium spectrum,

$$\frac{F_\mathrm{A}(E)}{E}=c_A\frac{F_4(E)}{E}(\mathrm{A}=15,27,56).$$

This assumption does not contradict the world data for the CNO and Ne-S nuclear groups but works a bit worse for the iron group. Nevertheless, a more sophisticated model would be unpractical since the corresponding correction would affect the secondary lepton fluxes by a negligible margin. A normalization using the relevant data yields $`c_{15}=0.068`$, $`c_{27}=0.026`$, $`c_{56}=0.0054`$.

In this paper we do not consider the effects of solar modulation and use the BESS+JACEE spectrum in all calculations. Therefore the predicted muon and neutrino fluxes are to some extent the maximum ones possible within our 1D cascade model.

## 4 Nucleon-nucleus and nucleus-nucleus interactions

Calculations with the earlier version of CORT  were based on semiempirical models for inclusive nucleon and light meson production in collisions of nucleons with nuclei, proposed by Kimel’ and Mokhov (KM)  and by Serov and Sychev (SS)  (see also  for the most recent versions). The KM model is valid for projectile nucleon momenta above $`4`$ GeV/c and for the secondary nucleon, pion and kaon momenta above 450, 150 and 300 MeV/c, respectively. Outside these ranges (that is mainly within the region of resonance production of pions) the SS model was used. Comparison of the KM model with the data and with some other models is discussed in .

Both models are in essence comprehensive parametrizations of the relevant accelerator data. In our opinion, the combined “KM+SS” model provides a rather safe and model-independent basis to the low-energy atmospheric muon and neutrino calculations. However it is not free of uncertainties. For the present calculation, the fitting parameters of the KM model for meson and nucleon production off different nuclear targets were updated using accelerator data not available for the original analysis . The values of the parameters were extrapolated to the air nuclei (N, O, Ar, C).444The old version of CORT used the renormalized $`N`$Be inclusive cross sections. The overall correction is less than 10-15% within the kinematic regions significant to atmospheric cascade calculations. Besides, energy-dependent correction factors were introduced into the model to tune up the output $`\pi ^+/\pi ^{}`$ ratio taking into account the relevant new data. This correction only slightly affects the total meson multiplicities and modifies the pion charge ratio within 5-9%.

The processes of meson regeneration and charge exchange ($`\pi ^\pm +\mathrm{Air}\pi ^{\pm ()}+X`$ etc.) are not of critical importance for production of leptons with energies of our interest and can be considered in a simplified way. Here we use a proper renormalization of the meson interaction lengths, which was deduced from the results of ref.  obtained for high-energy cascades.

The next important ingredient of any cascade calculations is a model for nucleus-nucleus collisions. The overwhelming majority of muon and neutrino flux predictions is based on the simplest superposition model (SM) in which the collision of a nucleus A, with total energy $`E_\mathrm{A}`$ and atomic weight $`A`$, against a target nucleus B, is treated as the superposition of $`A`$ independent collisions of nucleons with the target nucleus, each nucleon having an energy $`E_0=E_\mathrm{A}/A`$. This approximation is based on the hypothesis that, when the energy per nucleon of the projectile is much larger than the single nucleon binding energy, the $`A`$ nucleons interact incoherently. Generally speaking, this approach is rather far from reality and its applicability to the lepton flux calculations is not evident. A more accurate treatment of nucleus-nucleus collisions becomes especially important at low energies and for the regions or directions with high geomagnetic cutoffs. Indeed, the magnetic rigidity of a proton bounded in a nucleus is a factor $`A/Z`$ larger than that of a free proton of the same energy and therefore (considering the steepness of the primary spectrum) mainly nuclei are responsible for production of the low-energy secondaries.

Here we consider a modest generalization of a simple “Glauber-like” model used in . Namely, we write the inclusive spectrum of secondary particles $`c`$ ($`c=p`$, $`n`$, $`\pi ^\pm `$, $`K^\pm `$, $`K^0,\mathrm{}`$) produced in AB collisions as

$`{\displaystyle \frac{N_{\mathrm{AB}cX}}{x}}=`$ $`\xi _{\mathrm{AB}}^c\left[Z{\displaystyle \frac{N_{p\mathrm{B}cX}}{x}}+(AZ){\displaystyle \frac{N_{n\mathrm{B}cX}}{x}}\right]`$
$`+\left(1\xi _{\mathrm{AB}}^c\right)\left[Z\delta _{cp}+(AZ)\delta _{cn}\right]\delta (1x).`$ (3)
Here $`N_{N\mathrm{B}cX}/x`$ is the spectrum of particles $`c`$ produced in $`N\mathrm{B}`$ collisions ($`N=p,n`$) and $`\xi _{\mathrm{AB}}^c`$ is the average fraction of inelastically interacting nucleons of the projectile nucleus A. The term proportional to delta function describes the contribution of “spectator” nucleons from the projectile nucleus.

In the standard Glauber-Gribov multiple scattering theory the quantity $`\xi _{\mathrm{AB}}^c`$ is certainly independent of the type of inclusive particle $`c`$. On the other hand, it depends of the type of nucleus collision. Indeed, essentially all nucleons participate in the central AB collisions ($`\xi _{\mathrm{AB}}^c1`$)555Here we suppose for simplicity that the atomic weight of the projectile nucleus is not much larger than that of the target nucleus. while, according to the well-known Bialas–Bleszyński–Czyż (BBC) relation ,

$$\xi _{\mathrm{AB}}^c=\sigma _{N\mathrm{B}}^{\mathrm{inel}}/\sigma _{\mathrm{AB}}^{\mathrm{inel}}$$
(4)
for the minimum bias collisions.

To use (4) in a cascade calculation one should take into account that nucleons and mesons are effectively produced in nuclear collisions of different kind. Namely, the contribution from central collisions is almost inessential for the nucleon component of the cascade but quite important for light meson production. Thus one can expect that effectively $`\xi _{\mathrm{AB}}^{\pi ,K}>\xi _{\mathrm{AB}}^{p,n}`$.

From now on we use relation (4) for nucleon production by any nucleus while for meson production we put $`\xi _{\mathrm{He}\mathrm{Air}}^{\pi ,K}=\xi `$ and vary the parameter $`\xi `$ from 0.517 (the BBC value) to 0.710 (an experimental upper limit derived from the data on interactions of $`\alpha `$ particles with light nuclei ). Similar (within $`\pm 15`$%) variations of the parameters $`\xi _{\mathrm{A}\mathrm{Air}}^{\pi ,K}`$ for other nuclear groups lead to very small effects on the muon and neutrino fluxes; we do not discuss these effects below.

It is evident that the SM corresponds to the formal limit $`\xi _{\mathrm{AB}}^c=1`$ in eq. (4), if one simultaneously puts $`\sigma _{\mathrm{AB}}^{\mathrm{inel}}=\sigma _{N\mathrm{B}}^{\mathrm{inel}}`$ in the transport equations for the nuclei and mesons. Below, we use this approximation making efforts to represent the results of the Bartol code  which is based on the TARGET model for meson production and the SM for collisions of nuclei.

## 5 Numerical results: muon fluxes

In order to test the validity of our description of hadronic interactions and shower development, we compare in this section some numerical results for muon fluxes with data from several recent balloon-borne experiments. Moreover, in order to evaluate quantitatively theoretical uncertainties, we compare the calculations obtained by using different codes and/or interaction models. Namely, we consider numerical results obtained with:
i) the CORT code, using our “standard” (KM+SS) interaction model;
ii) the CORT code, using the TARGET model for $`\pi /K`$ meson production666Due to technical reasons it is difficult to adopt also the TARGET model for nucleon production; we thus describe the nucleon production by using the KM model everywhere. (below, we refer to this calculation as “CORT+TARGET”).
iii) the original Bartol Monte Carlo code  with the same (BESS+JACEE) fit for the primary spectra (“TARGET”).

In fig. 2 we compare the $`\mu ^{}`$ fluxes calculated with the three mentioned models for twelve atmospheric ranges with the data of the balloon-borne experiment CAPRICE 94 . The measurements were performed during ground runs in Lynn Lake ($`h=10^3`$ g/cm2) and during the continuous ascension of the balloon to its float altitude ($`h4`$ g/cm2). The nominal geomagnetic cutoff rigidity $`R_c`$ was about 0.5 GV and the detection cone opening was of about $`20^{}`$ around the vertical direction with the typical zenith angle of $`9^{}`$. The values of depths indicated in fig. 2 are the Flux-weighted Average Depths (FAD)  for eleven atmospheric ranges $`\mathrm{\Delta }h_i`$ ($`h_i<10^3`$ g/cm2).777There are also small (less than 3%) spreads of the FADs within some of the ranges $`\mathrm{\Delta }h_i`$, but these are completely negligible.

Figure 2 contains many data. In order to make it easily readable, let us explain the notation in detail. The solid curves represent the results of our calculations for each FAD. They are obtained with the standard CORT and with $`\xi `$ varying between the limits mentioned in section 4, that is within $`\pm 15`$%. The thickness of the curves reflects the indetermination in the parameter $`\xi `$. For obvious reasons, the muon flux uncertainty due to this indetermination is maximal at the top of the atmosphere (it is about 15% at $`h=3.9`$ g/cm2 and $`p=100`$ MeV/c) and becomes almost negligible for $`h100`$ g/cm2. The wide filled areas display the variations of the muon fluxes inside the ranges $`\mathrm{\Delta }h_i`$. It is important to note that the thickness of the bands is relatively small just for the region of effective muon and neutrino production that is in the neighborhood of the broad maxima of the muon flux ($`100300`$ g/cm2). This means that, in this region, the evaluation of the FAD cannot introduce any relevant uncertainty. Outside the region of effective production of leptons, the amplitude of the muon flux variations increases with decreasing muon momenta on account for the strong dependence of the meson production rate upon the depth at $`h100`$ g/cm2 and the growing role of the muon energy loss and decay at $`h300`$ g/cm2. The dashed curves show the results obtained using the original Bartol code and the dotted curves are the results of the CORT+TARGET model (we shall comment them later).

Figure 3 shows a comparison of the calculated differential momentum spectra of $`\mu ^+`$ and $`\mu ^{}`$ with the ground level data obtained in the experiments CAPRICE 94  ($`h=1000`$ g/cm2, $`R_c=0.5`$ GV) and CAPRICE 97  ($`h=886`$ g/cm2, $`R_c=4.2`$ GV). The detection cone in both experiments was the same as described above. The calculation are done using the KM+SS interaction model. Variation of $`\xi `$ leads to a small ($`3`$%) effect. One sees quite a good agreement between the theory and data.

For comparison, the result obtained using the SM is also shown. The data below $`10`$ GeV/c are precise enough in order to conclude that they disfavor the superposition model. Let us note that the calculated sea-level ($`h=10^3`$ g/cm2) $`\mu ^++\mu ^{}`$ flux is systematically lower than that obtained using the old version of CORT (see ref. ): the difference is about 18% for $`p=10`$ GeV/c and vanishes only at $`p80`$ GeV/c. The main reason of this difference is that the BESS+JACEE primary spectrum is essentially lower than the spectrum used in ref.  and in the earlier calculations  for energies below $`200`$ GeV/nucleon. Another reasons pertain to the changes in the inclusive cross sections and improved description of muon propagation.

In figs. 4 and 5 we present the muon charge ratio $`\mu ^+/\mu ^{}`$ as a function of muon momentum for different atmospheric depths and geomagnetic cutoffs. The data are from several balloon experiments: MASS 89 , MASS 91 , IMAX 92 , HEAT 94 , HEAT 95 (two data handlings of the same data sample given in  and ), CAPRICE 94  and CAPRICE 97 . The world average best fit, recently obtained in ref.  by an analysis of the data from many ground level experiments at $`p10`$ GeV/c, is also shown. The muon charge ratio, being sensitive to the primary chemical composition and to the $`\pi ^+/\pi ^{}`$ ratio at production, is a useful tool for testing the predicted $`\overline{\nu }/\nu `$ ratio.

Several important conclusions can be deduced by comparing the various data presented in figs. 2, 3, 4 and 5.

* There is a substantial agreement between the “standard” CORT predictions and the current muon data within wide ranges of muon momenta and atmospheric depths. In particular, the agreement is good for the region of effective production of leptons, in which the spread of the data (partially related to the procedure of determination of the FAD and mean muon momenta) is minimal. This provides an evidence for the validity of our model.
* The comparison between CORT and CORT+TARGET allows to elucidate the “hadronic” uncertainties in muon fluxes. The only difference between the two calculations is indeed the assumed meson production model and the model for nucleus-nucleus interactions. The $`\mu ^{}`$ fluxes obtained using the TARGET model systematically exceed the results obtained with our standard model. In order to have a quantitative estimate of the hadronic uncertainties, we note that the muon flux at $`p=1`$ GeV/c and $`h=(100300)`$ g/cm2 (i.e. the range where muons and neutrinos are effectively produced) increases by about 50% when the TARGET model is applied. For comparison, we remind that the accuracy of the CAPRICE 94 experiment in this region is about 15%.
* We can see that, except for very small depths, there is a serious disagreement between the CORT+TARGET and the original Bartol results. The disagreement grows fast with decreasing muon momentum and with increasing depth, and it is comparable (or even larger) to the hadronic uncertainties. Since the two calculations are obtained by using the same primary spectrum and the same meson production model, the observed difference can only be due to different features in the nuclear cascade development and muon propagation. On the other hand, the observed agreement at floating altitude, where such features are not important, assures the correctness of our use of the TARGET meson production model. Possible sources for the observed difference could be nucleon elasticity distributions, meson regeneration processes, muon energy losses and some more tiny features of the cascade models under consideration.

## 6 Numerical results: neutrino fluxes

In this section we present our numerical results for atmospheric neutrino (AN) flux. Due to geomagnetic effects, the AN spectra and angular distributions are quite different for different underground neutrino experiments. However, for the aims of present study it is enough to consider only one representative case. Let us confine ourselves to results pertinent to the Kamioka site.

Figure 6 shows the $`\nu _e`$, $`\overline{\nu }_e`$, $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ energy spectra averaged over all zenith and azimuth angles. The shaded areas are the results obtained with CORT using our standard interaction model. The widths of the areas indicate the uncertainty due to variations of the $`\xi `$ parameter. One sees that this uncertainty is at most 6% and thus it is negligible. The dashed curves correspond to the CORT+TARGET model, while the circles show the results by a new (preliminary as yet) 3D calculation based on the FLUKA Monte Carlo simulation package  (see also ). Note that the primary spectrum used in the latter calculation is also a parametrization of the BESS 98 data, which is however not identical to our BESS+JACEE fit. In fig. 7 we compare the zenith-angle distributions of $`\nu _e`$, $`\overline{\nu }_e`$, $`\nu _\mu `$ and $`\overline{\nu }_\mu `$, calculated with CORT (adopting $`\xi =0.685`$), CORT+TARGET and FLUKA; the distributions are averaged over azimuth angles and over seven energy ranges. The comparison allows to “highlight” the 3D effects which are highly dependent on neutrino energy and direction of arrival.

As a first point, we note that our standard model leads to the lowest AN flux at low energies. Our flux is also essentially lower than that predicted by Bartol group  and by Honda et al. 888Let us remind that the results of refs.  were obtained with different (higher) fluxes of the primary cosmic rays. and that are used in many analyses of the data from IMB, Fréjus, NUSEX, Kamiokande, Super-Kamiokande, SOUDAN 2 and MACRO.

We do not enter here in the complex problem of the interpretation of the AN anomaly in the light of our result. We want to remark, however, that, in the context of our model, it is difficult to increase the AN flux without spoiling the agreement with the current data on hadronic interactions, primary spectrum and muon fluxes. Let us emphasize that our low AN flux, when it is applied to the analysis of the underground neutrino experiments, results in some electron excess together with (or rather than) the muon deficit in the neutrino induced events.

Let us now discuss in detail the differences between our results and those obtained within CORT+TARGET and FLUKA. The quantitative comparison is given in tables 1 and 2. Table 1 shows the ratios of $`4\pi `$ averaged AN fluxes, obtained with CORT+TARGET and FLUKA to those with CORT for several energy ranges. In table 2 we tabulate the “flavor ratio”, defined by

$$R_\nu =\left(\nu _e+\frac{1}{3}\overline{\nu }_e\right)/\left(\nu _\mu +\frac{1}{3}\overline{\nu }_\mu \right),$$

where $`\nu _e`$, $`\overline{\nu }_e`$, etc. stand for the $`4\pi `$ averaged fluxes. This quantity is representative for the ratio of $`e`$ like to $`\mu `$ like single-ring contained events measured in water Cherenkov detectors and to “showers-to-tracks” ratio measured in iron detectors.

The comparison between CORT and CORT+TARGET predictions provides a direct determination of the “hadronic” uncertainties in the AN fluxes. We remind indeed that these two results are obtained by using the same primary fluxes and the same code; thus the observed differences is only due to the assumed meson production model. At low energies, the CORT+TARGET predictions drastically exceed those from CORT and even at energies about (2-2.5) GeV one can observe a sizeable difference between the two results. At $`E_\nu 0.5`$ GeV (tab. 1), the discrepancy is of the order of 40-50%. This shows that the hadronic uncertainties are quite relevant in the energy range which is responsible for the sub-GeV event rate in the Super-Kamiokande detector.

Figs. 6 and 7 show that CORT and FLUKA neutrino fluxes are in good agreement for any zenith angle at energies above $`1`$ GeV. On the other hand, CORT is systematically lower than FLUKA for $`E_\nu 1`$ GeV. At the energy $`E_\nu =0.5`$ GeV the FLUKA $`4\pi `$ averaged $`\nu _e`$, $`\overline{\nu }_e`$, $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ fluxes exceed the corresponding CORT results by about 15%, 27%, 18% and 18%, respectively. This discrepancy is only partially due to 3D effects (fig.7), which can account for an increase $`10\%`$, and it is most probably related to differences between the FLUKA and CORT hadronic interaction models.

Finally, as it is seen from table 2, the flavor ratio is essentially the same in all three models, except for very low energies. As one can conclude from fig. 7, the same is also true for the “up-to-down” asymmetry, that is for the ratio of upward-going ($`\mathrm{cos}\theta <0`$) and downward-going ($`\mathrm{cos}\theta >0`$) neutrino fluxes.

## 7 Conclusions

Let us summarize here the main points of our discussion.

* Calculations with a new version of the 1D code CORT, using the updated KM+SS interaction model and the BESS+JACEE primary spectrum are in agreement with the recent data on muon momentum spectra and charge ratios measured at different atmospheric depths and geomagnetic locations.
* The uncertainty due to indetermination of the parameters of our model for nucleus-nucleus interactions is comparatively small everywhere, except for high altitudes. Some experimental uncertainties (like those due to determination of the flux-weighted average depths) are negligible for the regions of effective generation of leptons and at ground level. Therefore the agreement between the predictions and the muon data obtained at these depths provides a conclusive evidence for validity of our approach.
* Our atmospheric neutrino flux is systematically lower than those used in the current analyses of the data from underground experiments, while the neutrino flavor ratio and “up-to-down” asymmetry are essentially the same as in all recent calculations. We stress that, in the context of our model, it is difficult to increase the AN flux without spoiling the agreement with the current data on hadronic interactions, primary spectrum and muon fluxes.

In order to estimate theoretical uncertainties in the muon and neutrino fluxes we compared our results with those obtained by using different codes and/or interaction models, with the following conclusions.

* By implementing two different interaction models within the same code, we have shown that the “hadronic” uncertainties are quite large at low energies for both muon and neutrino fluxes. Specifically, by using the TARGET meson production model instead of the KM+SS, we obtained a sizeable increase of neutrino fluxes at energies up to $`E_\nu 3`$ GeV. At the energy $`E_\nu 0.5`$ GeV representative for the sub-GeV events in the Super-Kamiokande detector, one observes a 40-50% increase of the AN fluxes.
* By comparing two different codes implementing the same interaction model (TARGET), we have shown that other sources of theoretical uncertainties are also relevant. Namely, a sizeable difference exists between the muon fluxes predicted by the original Bartol code and those obtained by implementing the TARGET interaction model within the CORT code. We believe that the observed disagreement is due to the differences in the calculation of nuclear cascade, muon energy loss and decay.
* A comparison of the $`4\pi `$ averaged AN fluxes calculated with CORT and FLUKA shows rather good agreement at $`E_\nu >0.70.8`$ GeV (the difference is typically less than 15%) and a sizable disagreement at lower energies. A comparison of neutrino zenith-angle distributions, predicted by CORT and FLUKA, suggests that the main sources of the disagreement at low energies are the differences in the adopted interaction models and 3D effects.

Finally we conclude that the AN spectra and angular distributions predicted with CORT using the BESS+JACEE primary spectrum and KM-SS interaction model are quite safe at energies above $`1`$ GeV while the low-energy range (which is of special importance, as a source of background to the proton decay experiments) is still uncertain.

We wish to acknowledge helpful discussions with G. Battistoni, M. Boezio, M. Circella and S. Coutu. We also thank them for providing us with their unpublished results. We are grateful to T. Gaisser and T. Stanev for making available to us their TARGET Monte Carlo package.


text:
# UTEXAS-HEP-01-20 MSUHEP-03161 Ultrahigh-Energy Neutrino-Nucleon Cross Sections and Perturbative Unitarity

## I Introduction

The phenomenology of ultrahigh-energy (UHE) neutrinos and their detection  depends on the neutrino-nucleon total cross section, which has been calculated in the standard model in Refs. . A striking feature of all of these predictions is the continued power-law-like growth of the cross sections with $`E_\nu `$ at the highest energies. This rise with $`E_\nu `$ is directly related to the very-low-$`x`$ behaviour of the nucleon parton distribution functions (PDF’s), which must be extrapolated below the regime of the current HERA data . It has even been argued  that cosmic ray events will soon require a cross section even larger than that given by the standard calculations.

On the other hand, unitarity relates the total scattering cross section to the forward scattering amplitude. For neutrino-nucleon scattering, in contrast to electron-nucleon scattering, the forward amplitude can be determined and used to bound the total cross section. This bound, based on lowest order perturbation theory in the weak coupling, is independent of the PDF’s and implies that the inelastic cross section cannot rise indefinitely. In fact, it cannot rise for more than about two decades in energy above the region now covered by HERA. Thus, either there will need to be a dramatic change in the way parton distribution functions scale not too far above the HERA data, or the assumption of weak perturbation theory is wrong. In either case the inelastic cross section will likely be very different at high $`E_\nu `$ than what is predicted in the references above.

In Section 2 we discuss the relation between the cross section at high $`E_\nu `$ and the low $`x`$ behavior of the distribution functions. In Section 3 we derive our bound and Section 4 contains conclusions.

## II Ultrahigh-Energy Neutrino-Nucleon Cross Sections

At ultrahigh neutrino energies the total cross section is completely dominated by its deep inelastic component. The deep inelastic scattering (DIS) cross section can be evaluated within the QCD-improved parton model, employing the universal parton distribution functions (PDF’s) of the nucleon  which are derived from fitting photon-exchange dominated DIS measurements. Most important for high energies are the recent HERA low-$`x`$ measurements , because UHE neutrino cross sections receive a dominant contribution from the ultrasmall-$`x`$ region of DIS – even below the range covered by HERA.

The total neutrino-nucleon cross section is given by

$$\sigma _{\mathrm{tot}}^{\nu NX}(s=2M_NE_\nu )=_0^1𝑑x_0^{sx}𝑑Q^2d^2\sigma ^{\nu NX}/𝑑x𝑑Q^2,$$
(1)
with the standard (weak) DIS cross section of the form Sub-leading QCD corrections are implicitly included in the evolution of $`q(x)=q(x,Q^2),\overline{q}(x,Q^2)`$; explicit $`𝒪(\alpha _s^1)`$ corrections to Eq. (1) are negligible $`𝒪(12\%)`$.

$$\frac{d^2\sigma ^{\nu NX}}{dxdQ^2}=\frac{G_F^2}{\pi }\left(\frac{M_{W,Z}^2}{Q^2+M_{W,Z}^2}\right)^2\left[q(x,Q^2)+\overline{q}(x,Q^2)(1y)^2\right],$$
(2)
where $`y=Q^2/(xs)`$. The effective quark and anti-quark distributions $`q(x,Q^2)`$ and $`\overline{q}(x,Q^2)`$ include the appropriate (electro-)weak couplings for charged or neutral current DIS. In what follows we will always assume an isoscalar nucleus $`[N=(p+n)/2]`$. Note, that Eq. (1) is well behaved only for neutrino DIS where – contrary to the case of photon exchange – the heavy vector meson propagator in Eq. (2) is non-singular  i.e. imposing a “$`Q>`$few GeV” cut-off on the integral in (1) changes its value only marginally. at $`Q^2=xys=0`$. Accordingly, the following arguments are tied to the weak interaction framework, with no analogue for the more familiar photon exchange process. For explicit calculations we restrict ourselves to neutrino, rather than anti-neutrino, scattering, with the understanding that sea-quark dominance in the UHE limit does not discriminate between neutrinos and anti-neutrinos.

Recent evaluations of UHE cross sections have extrapolated PDF’s below $`x<10^4`$ using power laws , radiative renormalization group behaviour  or BFKL-type resummations . Within the accuracy required for astrophysical phenomenology these approaches agree and predict a steep rise of the low-$`x`$ parton density functions. Any flattening of the PDF’s at low $`x`$ due to gluon recombination is naively expected to be absent at the UHE scale $`QM_W`$ due to the higher twist nature $`𝒪(1/Q^2)`$ of such recombination channels .

To gain a feeling for the interplay between the UHE cross section and the low-$`x`$ behaviour of the PDF’s we can perform the $`Q^2`$ integral under the approximation that we neglect the evolution of the PDF’s. We find\**\**\**We have verified that the effects of evolution induce corrections of only about 20%.

$$\sigma _{\mathrm{tot}}^{\nu NX}(s)\frac{M_{W,Z}^2G_F^2}{\pi }_0^1𝑑x\left[q(x)\left(\frac{\widehat{s}}{1+\widehat{s}}\right)+\overline{q}(x)\left(\frac{2}{\widehat{s}}+12\left(\frac{1+\widehat{s}}{\widehat{s}^2}\right)\mathrm{ln}(1+\widehat{s})\right)\right],$$
(3)
where $`\widehat{s}=sx/M_{W,Z}^2`$ and the PDF’s are evaluated at some fixed scale (i.e. $`q(x)=q(x,M_{W,Z}^2)`$, $`\overline{q}(x)=\overline{q}(x,M_{W,Z}^2)`$). From (3) it is straightforward to show that assumed behaviour of the PDF’s at $`x0`$ determines the $`E_\nu \mathrm{}`$ asymptotic behaviour of the total cross section : A power-type $`x0`$ behaviour $`xq(x)x^\beta `$  of the sea-quark PDF translates into a power-type rise with energy $`\sigma _{\mathrm{tot}}^{\nu NX}E_\nu ^\beta `$. Renormalization group evolution from a flat low-scale input  gives rise to a slower, but still significant, growth with energy. This is easiest to see in the double-leading-logarithm approximation , in which the cross section can be seen analytically  to rise less fast than any power but faster than a logarithm. Although each of the different approximations of the small-$`x`$ behaviour gives different quantitative predictions of the $`E_\nu \mathrm{}`$ asymptotics, we observe in Fig. 1 that they agree to a good approximation up to energies as high as $`E_\nu =10^{12}\mathrm{GeV}`$. In particular, they all rise continuously with energy.

This rise, which is a direct consequence of the $`x0`$ behaviour of $`q(x),\overline{q}(x)`$, reflects mainly the HERA measurements of $`F_2^\gamma xe_q^2(q+\overline{q})(x)`$ at small-$`x`$. The question is, therefore, up to what energy is the rise of $`\sigma _{\mathrm{tot}}`$ strictly required by these data? The $`ep`$ collider HERA runs at a center of mass energy of $`\sqrt{s}300\mathrm{GeV}`$ corresponding to an equivalent neutrino energy of $`E_\nu ^{\mathrm{HERA}}5\times 10^4\mathrm{GeV}`$ . However, leaving theoretical extrapolations aside, $`\sigma _{\mathrm{tot}}`$ for $`E_\nu >E_\nu ^{\mathrm{HERA}}`$ is not completely unknown because HERA probes the differential $`d\sigma /dxdQ^2`$ in Eq. (2) down to much lower $`x`$ (corresponding to higher $`E_\nu `$) than the average $`\mathrm{log}_{10}x1.5`$ in the integrated  $`\sigma _{\mathrm{tot}}(E_\nu ^{\mathrm{HERA}})`$, albeit at a lower scale $`Q<M_{W,Z}`$ for the $`\gamma ^{}`$-exchange process. Hence, the maximal $`E_\nu `$ up to which the evaluation of $`\sigma _{\mathrm{tot}}`$ according to Eqs. (1) and (2) is fixed by HERA depends on the extent to which we trust the NLO QCD scale-evolution of DIS structure functions. Recent estimates  of NNLO corrections to parton evolution confirm a trustworthy stability of the full singlet evolution at least down to $`x10^4`$ for $`Q^210\mathrm{GeV}^2`$. This latter kinematical range is well covered by HERA and the corresponding data are incorporated in the PDF sets of Refs. . Assuming nothing more††††††I.e. assumptions which are intrinsic to the perturbative evaluation in Eq. (1). than the universality of the PDF’s and the validity of the RGE within $`10^4x,\mathrm{\hspace{0.17em}10}\mathrm{GeV}^2<Q^2<M_W^2`$, we can then consider $`\sigma _{\mathrm{tot}}`$ – if evaluated according to Eqs. (1) and (2) – as effectively covered by HERA data up to neutrino energies where $`\sigma _{\mathrm{tot}}`$ becomes sensitive to $`x<10^4`$. This sets in‡‡‡‡‡‡For illustration, see Fig. 3 in  or Fig. 3 in . smoothly above $`E_\nu 10^6\mathrm{GeV}`$, which we indicate by the vertical line in Fig. 1. Note, this is a very conservative estimate of the impact of the HERA data. Were we to trust the PDF-mediated mapping of the HERA measurements onto the neutrino cross-section down to $`x3\times 10^5`$, and rely on the evolution from $`Q^2=1.5\mathrm{GeV}^2`$ to $`M_W^2`$, then the most recent HERA data  would cover $`\sigma _{\mathrm{tot}}`$ at $`E_\nu 10^8\mathrm{GeV}`$ except for a correction of $`25\%`$ from $`x<3\times 10^5`$. As we will deduce in the next section, this would suggest HERA $`\gamma `$-exchange data are already probing a perturbative unitarity bound in UHE neutrino scattering. However, for now, we prefer not to speculate about shifting the vertical line in Fig. 1 to the right. A detailed statistical analysis, though, might do better and provide stringent correlations between HERA data and $`\sigma _{\mathrm{tot}}`$ at even higher $`E_\nu `$. Note, our limitation $`x>10^4`$ is also a safe condition for the absence of recombination corrections  at HERA scales which are, accordingly, even more negligible at the UHE scale $`Q=M_{W,Z}`$.

## III Perturbative Unitarity

Using unitarity of the $`S`$-matrix, we can relate the total $`\nu N`$ cross section in (1) to the imaginary part of the neutrino-nucleon forward scattering amplitude:

$$\sqrt{\lambda }\sigma _{\mathrm{tot}}^{\nu NX}(s)=\mathrm{Im}\left[T_{\{\nu N,\nu N\}}(s,t=0)\right],$$
(4)
where $`\lambda =(sM_N^2)^2`$ and $`s,t`$ are the standard Mandelstam variables. The elastic amplitude $`T`$ is related to the elastic cross section

$$\frac{d\sigma _{\mathrm{el}}^{\nu N\nu N}}{dt}=\frac{1}{16\pi \lambda }\left|T_{\{\nu N,\nu N\}}(s,t)\right|^2.$$
(5)
Combining Eqs. (4) and (5) gives a general limit on $`\sigma _{\mathrm{tot}}`$ in Eq. (1)

$$\frac{d\sigma _{\mathrm{el}}}{dt}|_{t=0}=\frac{1}{16\pi }\left[\frac{\left(\mathrm{Re}[T_{\{\nu N,\nu N\}}(s,0)]\right)^2}{\lambda }+\sigma _{\mathrm{tot}}^2\right]\frac{1}{16\pi }\sigma _{\mathrm{tot}}^2.$$
(6)
We note that the inequality Eq. (6), as derived, holds strictly for each spin and isospin state of the nucleon; however, it is straightforward to show that\**\**\**The difference between the left and right hand sides of Eq. (7) is the variance of elementary probability theory and therefore positive.

$$\frac{1}{4}\sigma _{\mathrm{tot}}^2\left[\frac{1}{4}\sigma _{\mathrm{tot}}\right]^2,$$
(7)
where the sum is over spin and isospin of the nucleon. Thus, we can just as well use the spin/isospin averaged cross sections in Eq. (6)\*†\*†\*†In fact, due to constraints on spin and isospin asymmetries in DIS scattering, one would expect the difference between the right-hand and left-hand sides of Eq. (7) to be negligible in the high energy limit..

The inequality in Eq. (6) is a standard statement of the optical theorem and can be found in many textbooks. It follows strictly from the positivity of $`(\mathrm{Re}T)^2`$ and does not rely on any perturbative expansion. We now consider its implications in the context of a perturbative calculation in the weak coupling. Expanding the elastic amplitude to lowest power in the weak coupling $`g`$ and using the most general $`Z`$-nucleon coupling (see Fig. 2), we obtain

$`T`$ $`=`$ $`\left({\displaystyle \frac{g}{4\mathrm{cos}\theta _W}}\right)^2\overline{u}_\nu (k^{})\gamma ^\alpha (1\gamma _5)u_\nu (k){\displaystyle \frac{1}{tM_Z^2}}`$ (10)
$`\times \overline{u}_N(p^{})[\gamma _\alpha f_1(q^2)+i\sigma _{\alpha \beta }q^\beta f_2(q^2)+q_\alpha f_3(q^2)+\gamma _\alpha \gamma _5g_1(q^2)`$
$`+i\sigma _{\alpha \beta }\gamma _5q^\beta g_2(q^2)+q_\alpha \gamma _5g_3(q^2)]u_N(p)+𝒪(g^4).`$
In the forward direction only the form factors $`f_1(0)`$ and $`g_1(0)`$ contribute. Using the standard model $`Z`$-boson current and isospin symmetry, the values of the form factors at $`q^2=0`$ can be expressed in terms of the nucleon isospin operator $`T_3`$ and charge operator $`Q`$ as $`f_1(0)=(2T_34\mathrm{sin}^2\theta _WQ)`$ and $`g_1(0)=2T_3g_A`$, where $`g_A=1.27`$ is the measured value  of the axial vector coupling constant. For an isoscalar nucleus we obtain

$$\frac{d\sigma _{\mathrm{el}}^{(2)}}{dt}|_{t=0}=\frac{G_F^2}{8\pi }\left[\frac{1+(14\mathrm{sin}^2\theta _W)^2}{2}+g_A^2\right].$$
(11)
Accordingly, the inequality Eq. (6) gives

$$\sigma _{\mathrm{tot}}9.3\times 10^{33}\mathrm{cm}^2.$$
(12)
The value on the right-hand side of Eq. (12) is shown as a horizontal line in Fig. 1. We observe that the inequality is violated for $`E_\nu 2\times 10^8\mathrm{GeV}`$. Note that the precise value of the bound is rather insensitive to the values of $`f_1(0)`$ and $`g_1(0)`$ because it depends on the square root of the factor in the square brackets of Eq. (11).

HERA information on PDF’s determines $`\sigma _{\mathrm{tot}}`$ for neutrino energies as high as $`E_\nu 10^6\mathrm{GeV}`$. Perhaps large changes in the PDF’s will set in within the range of $`x`$ appropriate for $`E_\nu 10^610^8\mathrm{GeV}`$ such that the bound is respected for larger $`E_\nu `$. These corrections would presumably arise as a softening of the gluon and sea quark distributions at small $`x`$.

Alternatively, one can note that the bound Eq. (12) is obtained by comparing an $`𝒪(g^4)`$ expression for $`\sigma _{\mathrm{tot}}`$ with an $`𝒪(g^2)`$ expression for the forward elastic scattering amplitude. Since unitarity is an order-by-order statement within perturbation theory, to be strictly rigorous we must also include the absorptive $`𝒪(g^4)`$ contribution to the forward scattering amplitude. In this way, unitarity in the ultrahigh energy limit can in principle be restored, but at the cost of large $`𝒪(g^4)`$ corrections to the forward amplitude.

For illustration, we write the elastic amplitude symbolically as

$$T=g^2T^{(2)}+g^4T^{(4)}+\mathrm{}.$$
(13)
where we have calculated the $`T^{(2)}`$ part above and the $`T^{(4)}`$ term results from the exchange of a pair of virtual $`Z^0`$’s or $`W^\pm `$’s in the elastic scattering, as illustrated in Fig. 3. Using this notation

$$\frac{d\sigma _{\mathrm{el}}}{dt}|_{t=0}g^4T_{}^{(2)}{}_{}{}^{2}+2g^6T^{(2)}\mathrm{Re}[T^{(4)}]+g^8\left[\left(\mathrm{Re}[T^{(4)}]\right)^2+\left(\mathrm{Im}[T^{(4)}]\right)^2+\mathrm{}\right]+\mathrm{},$$
(14)
and

$$\left(\sigma _{\mathrm{tot}}\right)^2g^8\left(\mathrm{Im}[T^{(4)}]\right)^2+\mathrm{}.$$
(15)
From these equations we see that the unitarity condition Eq. (6) can be satisfied; however it requires that the higher order $`g^4T^{(4)}`$ term become larger than the lower order $`g^2T^{(2)}`$ term in the high energy limit. So, in this approach, calculating to lowest order in the weak coupling $`g`$ is not a good approximation for the elastic differential cross section at high energy. But if the perturbation expansion is unreliable in Eq. (14), then it is also reasonable to expect higher order corrections to Im$`T`$ to be larger than the term shown, both in Eq. (14) and in Eq. (15). Since this is the term on which all calculations have been based it implies that these calculations may not be reliable at very high energies. If the perturbation expansion doesn’t work in Eq. (14) then we must question its validity in Eq. (15).

## IV Conclusions

We compared perturbative predictions for neutrino-nucleon total cross sections at ultrahigh neutrino energies to a unitarity bound derived from the corresponding elastic neutrino scattering amplitude in the forward direction. At the non-perturbative level the bound is absolutely stringent. A lowest order expansion in powers of the weak coupling leads to a bound which is saturated at an energy surprisingly close to what is effectively covered by HERA measurements. Thus, for the largest energies relevant to neutrino astronomy, sizable PDF corrections could reside in $`\sigma _{\mathrm{tot}}`$. Alternatively, the large apparent $`𝒪(g^4)`$ total cross section should manifest itself through the corresponding term in the elastic amplitude as well. But if the perturbation expansion is not valid for the elastic amplitude how can we trust it for the inelastic cross section? In this case also $`\sigma _{\mathrm{tot}}`$ could be very different for large $`E_\nu `$ than the curves shown in Fig. 1.

We have ignored the possibility that unitarity is satisfied by a new interaction – physics beyond the standard model . But even in this alternative our basic conclusions remain the same: something dramatic must happen at energies close to the current energy and the existing calculations of $`\sigma _{\mathrm{tot}}`$ should not be trusted at large neutrino energy.

If the bound is correct, and the PDF’s change at low $`x`$, then $`\sigma _{\mathrm{tot}}`$ will be smaller at high $`E_\nu `$ than what is shown in Fig. 1. If perturbation theory has broken down then $`\sigma _{\mathrm{tot}}`$ could be anything - perhaps even large enough to explain the cosmic ray data without the need for new physics. Or perhaps the values shown in Fig. 1 could turn out to be correct; maybe weak perturbation theory only breaks down for the elastic cross section. Without a better knowledge of the PDF’s at low $`x`$, there is no way to know.

###### Acknowledgements.

It is a pleasure to thank Karol Lang, Wu-Ki Tung and Scott Willenbrock for helpful conversations and G. Domokos, M. Glück, and S. Kovesi-Domokos for their comments. This research was supported in part by the National Science Foundation under Grant PHY-0070443 and by the United States Department of Energy under Contract No. DE-FG03-93ER40757.