text:
# LISA Sources in Globular Clusters

## 1 Introduction

Binary stars will be the most common continuous sources of gravitational radiation (GR) for LISA. The detection of these sources is enhanced for systems with short orbital periods and large component masses, favoring double compact stars, i.e., binaries containing white dwarf (WD), neutron star (NS), or black hole (BH) components. Globular clusters will be excellent hosts as they are known to contain large populations of compact objects in binaries formed through dynamical interactions (e.g., Bailyn 1995; Hut et al 1992). While some of these systems have been detected directly through their electromagnetic radiation (e.g., as radio pulsars, X-ray or UV sources), most are thought to be invisible in conventional astronomy, but they could be very bright in gravitational radiation. For example, a large, dense globular cluster such as 47 Tuc may contain as many as $`10^4`$ NS, of which $`10^3`$ may be active radio millisecond pulsars, but only $`20`$ are bright enough to have been detected in pulsar surveys (Camilo et al 2000). Most of the detected pulsars are in binaries with low-mass WD companions.

If LISA can detect and identify some of these sources within specific globular clusters, then gravitational wave astronomy will be able to contribute significantly to our understanding of binary star evolution and stellar dynamics in dense cluster systems. Many Galactic globular clusters are close enough that GR emission from double degenerate binaries may be detected by LISA with a signal-to-noise (S/N) ratio $`>2`$. At this level, LISA will have an angular resolution of a few square degrees. Individual clusters will be smaller than this, but the typical separation between clusters is much larger, so that LISA will be able to separate sources in different clusters. One of the primary sources of noise in the frequency band of cluster binaries is the Galactic disk population of close white-dwarf binaries, which is expected to produce a confusion-limited signal above the noise level of LISA (Hils and Bender 2000, Schneider et al2000). It may be possible to avoid this noise by looking specifically for signals originating out of the plane of the Galaxy, such as those arising from globular clusters. In this paper, we discuss the possible source populations and we present some simple methods for crudely estimating their numbers.

Since the angular resolution of LISA depends upon measuring the Doppler shift in frequency as a phase variation during the motion of the detector about the sun, the wavelength of the signal must be less than the size of LISA’s orbit. This places an upper limit on the binary orbital period, $`P_{\mathrm{orb}}=2/f_{\mathrm{gw}}=2\lambda /c4AU/c2000`$s. In addition, we require that the S/N ratio be $`>2`$ for the binary to be detectable with the required positional accuracy. By coincidence, this also imposes a maximum period $`P_{\mathrm{max}}2000`$s above which the confusion-limited noise from Galactic binaries rises significantly. Finally, we require that the systems be detached so that their period evolution is solely determined by GR, providing a “clean” signal. This should also improve the extraction of the signal from the confusion-limited background. These requirements limit the viable sources to double compact objects. We will not discuss other standard LISA sources such as contact binaries containing two main-sequence stars or a main-sequence star and a WD, as their signals will not be strong enough or clean enough. Integrating the in-spiral of a binary driven by GR alone from $`P_{\mathrm{max}}`$ down to contact, we can write the total length of time spent in the LISA band as

$$T_{\mathrm{LISA}}4\times 10^6\mathrm{yr}\left(\frac{M_{\mathrm{chirp}}}{M_{}}\right)^{5/3}\left(\frac{P_{\mathrm{max}}}{2000\mathrm{s}}\right)^{8/3},$$
(1)
where $`M_{\mathrm{chirp}}=(m_1m_2)^{3/5}/(m_1+m_2)^{1/5}`$ is the chirp mass, given in terms of the binary component masses $`m_1`$ and $`m_2`$.

## 2 Theoretical Estimates

Numerical simulations of globular cluster dynamics can be used to try to estimate the numbers and properties of different classes of binary systems from first principles. Here we illustrate this approach from three recent studies using different numerical methods.

### 2.1 $`N`$-body simulations

The dynamical simulations of Portegies Zwart and McMillan (2000, hereafter PZM) were developed to estimate BH-BH merger rates for LIGO and VIRGO. In their $`N`$-body simulations of globular cluster evolution, a component of $`10M_{}`$ primordial BH quickly undergoes mass segregation, forming an effectively isolated subsystem at the cluster core in $`10^9`$yr. Through dynamical interactions, hard BH binaries are then formed and ejected from the cluster and the entire BH subsystem evaporates typically within $`5\times 10^9`$yr. The basic conclusion of this work is that globular clusters today should no longer contain any of their primordial BH. This explains naturally why no BH X-ray binary has ever been detected in a globular cluster, even though clusters contain many bright NS X-ray binaries (see, e.g., Deutsch et al 2000). A large number of ejected BH binaries will coalesce within a Hubble time, and these may provide a significant event rate for ground-based interferometers.

However, the last BH binary remaining in the cluster core may not be ejected, since it can now only interact with the (much less massive) normal stars. The dynamical simulations of PZM suggest that the initial separations of typical BH binaries should be $`30300R_{}`$. The eccentricities are typically large, and roughly half of these binaries will coalesce through GR within the remaining $`t_{\mathrm{coal}}5\times 10^9`$yr of the cluster evolution to the present (it is easy to show that, for coalescing binaries, further hardening through dynamical interactions with normal cluster stars is negligible, and that ejection from the cluster is extremely unlikely). Unfortunately, with $`T_{\mathrm{LISA}}10^5`$yr for two $`10M_{}`$ BH (see eq. 1), even if we assume optimistically that $`100`$ globular clusters in our Galaxy retained a coalescing BH binary, we still conclude that the number of BH-BH systems in the LISA band today is $`100T_{\mathrm{LISA}}/t_{\mathrm{coal}}1`$. Thus BH-BH binaries in clusters are not promising sources for LISA.

### 2.2 Encounter Rates

Davies (1995) calculated the production rates of various binaries through dynamical interactions using a combination of scattering experiments and static models of globular clusters. His results can be used to estimate the numbers of NS-NS, NS-WD, and WD-WD binaries with $`P_{\mathrm{orb}}2000`$s (Benacquista 1999).

The original dynamical simulation calculates three-body interaction cross-sections for a variety of interactions and outcomes and then determines encounter rates for a number of model globular clusters. Production rates of binaries were inferred from the cross-sections applied to the encounters. As a side result, the merger rates through GR induced in-spiral of binaries were produced for 15 Gyr runs. The runs were also extended until all binaries had merged. This generally took another 15 Gyr. With these two results, a merger rate can be determined through

$$\eta =\frac{(n_{\mathrm{}}n_{15})}{15\mathrm{Gyr}}$$
(2)
where $`n_{\mathrm{}}`$ is the number of mergers after all binaries have merged and $`n_{15}`$ is the number which have merged after the first 15 Gyr. The evolution of the orbital period can be determined from $`\dot{P}=k_oP^{5/3}`$ with

$$k_o=\frac{96}{5}(2\pi )^{8/3}\frac{G^{5/3}}{c^5}M_{\mathrm{chirp}}^{5/3}.$$
(3)
Combining these results in a number distribution given by

$$dn=\frac{\eta }{k_o}P^{5/3}dP.$$
(4)
Using the numbers found in Davies (1995) and scaling up to the entire globular cluster system, with a total number of stars $`N_{\mathrm{gc}}10^{7.5}`$, we find 16 NS-NS binaries, 125 NS-WD binaries, and 176 WD-WD binaries in the LISA band.

### 2.3 Monte Carlo Simulations

Dynamical Monte Carlo simulations can be used to study self-consistently the evolution of binary populations within evolving globular cluster models (Rasio 2000). Rasio et al 2000 used this approach to study the formation and evolution of NS-WD binaries, motivated by the recent detection of a large number of millisecond radio pulsars with WD companions in 47 Tuc (Camilo et al 1999). The simulations include an evolving core of the background cluster, stellar evolution, mass segregation, and dynamical 3-body interactions. In the process of generating the appropriate populations of binary millisecond pulsars, the simulations also produce populations of NS-WD binaries in the period range of interest for LISA. There are two such binaries in the simulation of 47 Tuc (Rasio et al 2000). From this admittedly small sample it is possible to place an order of magnitude estimate on the number of such binaries in the globular clusters system. Assuming that the 2 NS-WD binaries (presently in the LISA band) for the $`10^6M_{}`$ in 47 Tuc is consistent throughout the globular cluster system, we would predict a total of $`60`$ such binaries for the $`10^{7.5}M_{}`$ in all clusters. Note that this number is within a factor of two of the number estimated above from encounter rates.

## 3 Semi-Empirical Estimates

Alternatively, we can estimate the numbers of the various possible LISA sources by trying to relate them to certain types of binaries already observed in electromagnetic radiation. Here the best candidates are NS-WD systems, also observable as X-ray binaries, and NS-NS systems, also observable as binary radio pulsars.

### 3.1 NS-WD Binaries

The small number of observed ultra-compact X-ray binaries (UXBs) in globular clusters (Deutsch et al 2000) can be used to estimate empirically the number of NS-WD binaries in the LISA band. UXBs are bright enough that we essentially see all of them in the Galactic globular cluster system. Assuming that they originate from the GR-driven in-spiral of NS-WD binaries all the way to contact (Rasio et al 2000) and that the process is ongoing, we can estimate the number of progenitor systems in the LISA band by looking at the time spent by a system in both phases. Let $`N_\mathrm{X}`$ be the number of UXBs, and $`T_\mathrm{X}`$ their typical lifetime. Also, let $`N_{\mathrm{LISA}}`$ be the number of progenitors in the LISA band, with the corresponding $`T_{\mathrm{LISA}}`$ given by eq. 1. Then, if the process is stationary, we must have

$$\frac{N_\mathrm{X}}{T_\mathrm{X}}\frac{N_{\mathrm{LISA}}}{T_{\mathrm{LISA}}}.$$
(5)
Of the 12 known bright X-ray binaries in clusters, 3 are known to have very short orbital periods ($`<1`$hr) and 3 are known to have long orbital periods (Deutsch et al1999). Thus, $`3N_\mathrm{X}9`$. The lifetime $`T_\mathrm{X}`$ is rather uncertain, depending on a number of poorly understood processes. A standard treatment of mass transfer driven by GR only would give $`T_\mathrm{X}10^9`$yr, but other effects such as irradiation or tidal heating of the companion may well shorten this to $`T_\mathrm{X}10^7`$yr, which would then also provide consistency with the birthrate of binary millisecond pulsars (Applegate and Shaham 1994; Tavani 1991; Rasio et al2000). For a typical system containing a $`1.4M_{}`$ NS and a $`0.5M_{}`$ WD companion, eq. 1 gives $`T_{\mathrm{LISA}}7\times 10^6`$yr. Adopting the most optimistic values of $`N_\mathrm{X}`$ and $`T_\mathrm{X}`$ we arrive at $`N_{\mathrm{LISA}}110`$, about an order of magnitude smaller than predicted by dynamical simulations.

### 3.2 NS-NS Binaries

One NS-NS binary (M15C) has been detected directly as a binary radio pulsar, out of a total of $`40`$ radio pulsars known in all Galactic globular clusters. A proper statistical analysis taking into account the various selection effects in pulsar surveys is well beyond the scope of this paper. However, since the dynamical processes that created M15C are thought to happen quite generically in all dense clusters (see, e.g., Phinney and Sigurdsson 1991), and since selection effects, if anything, tend to make pulsars in close binaries more difficult to detect than other types of radio pulsars (Johnston and Kulkarni 1991), it may not be unreasonable to estimate the total number of NS-NS binaries in globular clusters as $`N_{\mathrm{NS}}/40`$, where $`N_{\mathrm{NS}}`$ is the total number of NS in clusters. A strict lower limit on $`N_{\mathrm{NS}}`$ comes from estimates of the total number of active radio pulsars in clusters, giving $`N_{\mathrm{NS}}>10^4`$ (see, e.g., Kulkarni et al 1990). An upper limit comes from dynamical models of globular clusters suggesting that as much as $`5\%`$ of the total mass of a cluster may be in the form of $`1.4M_{}`$ dark remnants (some of which may be massive WDs, which, for our purposes, are indistinguishable from NS; see, e.g., Murphy et al 1998). Taking again the total mass in globular clusters $`10^{7.5}M_{}`$ (Prince et al 1991), this gives $`N_{\mathrm{NS}}<10^6`$. Adopting the coalescence time of M15C, $`t_{\mathrm{coal}}=2\times 10^8`$yr, as typical, and using $`T_{\mathrm{LISA}}=3\times 10^6`$yr from eq. 1 for two $`1.4M_{}`$ NS, we get a number $`1100`$ NS-NS binaries in the LISA band, in rough agreement with the predictions of dynamical simulations.

## 4 Detection Efficiency

From the above estimates, we see that there may be many compact binaries in the Galactic globular cluster system that will be detectable GR sources for LISA. In order to determine the efficiency of LISA for detecting these sources throughout the globular cluster system, we calculated the S/N ratio for a range of binary masses and frequencies from each globular cluster in the system. We begin by modeling the strain amplitude measured by LISA. A stationary detector will measure a strain amplitude found using the point-mass quadrupole formula for a circularized binary given by

$$h(t)=\frac{\sqrt{3}}{2}(A_+F^+\mathrm{cos}(2\pi ft)+A_\times F^\times \mathrm{sin}(2\pi ft)),$$
(6)
where $`F_+`$ and $`F_\times `$ represent the sensitivity of the detector to each polarization. The amplitudes are given by

$`A_+={\displaystyle \frac{G^{5/3}}{c^4d}}M_{\mathrm{chirp}}^{5/3}f^{2/3}(1+\mathrm{cos}^2i)`$ (7)
$`A_\times =2{\displaystyle \frac{G^{5/3}}{c^4d}}M_{\mathrm{chirp}}^{5/3}f^{2/3}\mathrm{cos}i.`$ (8)
The angle of inclination, $`i`$, is the angle that the line of sight makes with the orbital plane of the binary. Since LISA is in orbit about the sun and its orientation changes, the actual signal received by LISA is given by

$$h(t)=\frac{\sqrt{3}}{2}A(t)\mathrm{cos}\left[^t2\pi f(t^{})𝑑t^{}+\phi _p(t)+\phi _D(t)+\phi _0\right].$$
(9)
The Doppler shift in frequency due to LISA’s motion about the sun is incorporated in $`\phi _D(t)`$, and the changing orientation of LISA with respect to the source is incorporated in both $`\phi _p(t)`$ and $`A(t)=\sqrt{(A_+F^+)^2+(A_\times F^\times )^2}`$. The exact form of $`\phi _D`$, $`\phi _p`$, $`F^+`$, and $`F^\times `$ can be found in Cutler (1998). The S/N ratio $`\rho `$ is then found from

$$\rho ^2=4\underset{\alpha }{}_0^{\mathrm{}}\stackrel{~}{h}_\alpha ^{}(f)\stackrel{~}{h}_\alpha (f)\frac{df}{S_n(f)},$$
(10)
with the sum over $`\alpha `$ representing the sum over both semi-independent arms of LISA and the Fourier transform $`\stackrel{~}{h}`$ defined as

$$\stackrel{~}{h}_\alpha (f)=_{\mathrm{}}^{\mathrm{}}\mathrm{}^{2\pi \mathrm{}ft}h_\alpha (t)𝑑t.$$
(11)
Using globular cluster position and distance data from Harris (1996), we determined the S/N ratio for NS-WD binaries with WD masses in the range of 0.1 to 0.4 $`M_{}`$ and frequencies ranging from $`1`$mHz up to Roche lobe overflow, using the expression of Rappaport et al (1987), $`f_{\mathrm{max}}=43.3\mathrm{mHz}(M_{\mathrm{WD}}/M_{})`$. We note that the signal strength for a binary with $`i=0`$ is about twice that for one with $`i=\pi /2`$. The globular clusters were then classified by the strength of the detected signal for binaries with favorable or unfavorable angles of inclination. There are 7 clusters with some binaries having $`i=\pi /2`$ and $`\rho 5`$ and 17 for which only binaries with $`i=0`$ have $`\rho 5`$. The binaries associated with these clusters would have sufficient angular resolution to clearly identify the binary with its cluster. There are 72 clusters with $`i=\pi /2`$ and $`\rho 2`$ and 36 with only $`i=0`$ and $`\rho 2`$. These clusters could house binaries whose angular resolution would be a few square degrees, and the identification of the binary with its host cluster would be a little more suspect. There are only 15 clusters which are too distant to detect any binaries with $`\rho 2`$.

## 5 Conclusion

There is ample, direct and indirect, evidence for the existence of a population of compact binaries in globular clusters that will be detectable GR sources for LISA. The observed populations of X-ray binaries and millisecond radio pulsars clearly suggest a rich and interesting variety of dynamical processes that should produce large numbers of compact binaries. A variety of estimation techniques, both theoretical and semi-empirical, are in rough agreement to predict as many as $`100`$ NS-NS, NS-WD, and WD-WD binaries detectable by LISA in the Galactic globular cluster system. Although there are very large uncertainties in the exact numbers, the existence of these populations is not in doubt. LISA will be able to probe about 90% of Galactic globular clusters to search for binaries with orbital periods below 2000 s. Measurements of this population of binaries will provide valuable insight into the dynamical evolution of dense star clusters and the role that binaries play in them.

M.J.B. acknowledges support from NASA EPSCoR grant NCCW-0058 and Montana EPSCoR grant NCC5-240. S.P.Z. was supported by NASA through Hubble Fellowship grant HF-01112.01-98A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy. F.A.R. acknowledges support from NSF Grants AST-9618116 and PHY-0070918, and from NASA ATP Grant NAG5-8460.

## References


format: text

text:
# Dark matter distribution in the universe and ultra-high energy cosmic rays

## Introduction

Most of our universe is made of a dark elusive component, whose nature is still unknown. This mysterious component dominates the gravitational interactions in our world and determines the distribution of visible matter, that we associate with galaxies, clusters and superclusters of galaxies and their luminous constituents. We know from the measure of the values of cosmological parameters, that most of the dark matter must be of non baryonic nature. This important finding stimulated the blossoming of numerous proposals of new particles, in the context of theories of Grand Unification or extensions of the Standard Model. Some of the possibilities that have been proposed invoke the existence of super heavy particles, with mass comparable with the mass of the inflaton. These particles might be almost stable and contribute the whole of cold dark matter needed to explain cosmological and astrophysical observations. Independent of the nature of dark matter, astrophysical observations of the visible universe allow us to know how dark matter is distributed in space: primordial fluctuations in the density field result today in a structured distribution of the dark and visible components in the universe. Galaxies are part of this clustering, but similar local overdensities can be seen on larger scales, up to a few tens Mpc, where the universe starts being well described by a homogeneous density field.

An apparently unrelated mystery of modern physics is the existence of particles with energy in excess of $`10^{20}`$ eV in the cosmic radiation, the so-called ultra high energy cosmic rays (UHECRs). These particles are thought to be of extragalactic origin, bacause the galactic magnetic field is unable to confine particles with such energy. As pointed out in the pioneering papers by Greisen greis  and Zatsepin and Kuzmin kz , if the sources are distributed homogeneously in the universe, the process of photopion production off the photons of the cosmic microwave background (CMB) should cause a strong suppression in the spectrum of UHECRs above $`5\times 10^{19}`$ eV, what is now known as the Greisen-Zatsepin-Kuzmin (GZK) cutoff. Nevertheless we do detect particles with energy much in excess of this cutoff energy, and we have been unable to identify any straightforward astrophysical source in the direction of arrival of these events. Can the two mysteries, of dark matter and UHECRs be related? We discuss here two possible connections: in sections 2 and 3 we explore the consequences that a local clustering of the sources, driven by the corresponding dark matter clustering, have on the observed fluxes of UHECRs, in comparison with the case of homogeneous distribution. The initial prediction by Greisen and Zatsepin and Kuzmin, was in fact derived in the assumption of sources equally distributed in space. It was later shown by Berezinsky and Grigorieva beregri79  that a local overdensity by a factor $`10`$ over $`20`$ Mpc would make the problem of the existence of UHECRs less severe. The local density of sources can now be extracted from large catalogs, so that a realistic determination of the density field can finally be used in the calculation of the fluxes of UHECRs. We use here the PSCz and CfA catalogs for the purpose of calculating the density field. The fluxes of UHECRs are calculated numerically by montecarlo simulations.

In section 4 we consider a direct connection between dark matter and UHECRs, in the case that dark matter is composed of super heavy quasi-stable particles. Particles fulfilling these requirements have been recently discussed in kolb ; bere ; kuzmin ; kt  \[see also J. Ellis (these proceedings)\]. We will review here the current situation and stress some consequences of this model on the observable large scale and small scale anisotropies in the arrival directions of UHECRs.

## Measuring the Galaxy Density Field

In order to determine the effects of the inhomogeneity of the source distribution on the cosmic ray spectrum, we need to measure the galaxy density field, and in particular its dependence on redshift. In general, it may well be that the sources of UHECRs are of some type not directly related to other known objects like galaxies, but we assume here that the sources of UHECRs have a density field that is proportional to that of ordinary galaxies. We follow the approach presented in bbo00 . The galaxy density field is usually measured by selecting galaxies from an imaging sky survey and taking their redshifts. Almost invariably, the galaxies are selected to be brighter than some limiting flux $`f_{\mathrm{lim}}`$ in some band, expressed as an “apparent magnitude” $`m_{\mathrm{lim}}=2.5\mathrm{log}_{10}(f_{\mathrm{lim}}/f_0)`$, where $`f_0`$ is an arbitrary zero-point. For all (or for some random subsample) of the galaxies brighter than this, their spectra are taken and their redshifts $`z`$ are determined.

However, we cannot simply use the raw distribution of redshifts from such a flux-limited survey, regardless of the way the galaxies were selected. Here we describe the proper way to derive density fields from galaxy redshift surveys. We limit ourselves to measuring the density in redshift space, ignoring the effects of deviations from the Hubble law due to galaxy peculiar velocities. The subject of galaxy density fields dates back to Ref. davis82a . An educational recent review is that of Ref. strauss95a .

Our results will be based on two surveys. First, we consider the Center for Astrophysics Redshift Survey (CfA2; huchra95a ). This survey comprises about 10,000 galaxy redshifts, selected to be brighter than $`m=15.5`$, (approximately, a $`B`$-band magnitude). It covers an area which is about 17% of the whole sky. However, in order to evaluate the effects of the density field of galaxies on the cosmic ray spectrum, we really should probe the density field over nearly the whole sky. The best sample of galaxies to use for this purpose is the IRAS PSC$`z`$ Survey saunders00a , which consists of about 15,000 galaxies with infrared fluxes $`>0.6`$ Jy over about 84% of the sky.

A consequence of the flux limits in any survey is that at different redshifts, a different set of galaxy luminosities $`L`$ is observed, determined by the faintest luminosity observable at that redshift $`L_{\mathrm{min}}(z)`$. For an Euclidean metric, this luminosity is related to the flux limit by $`L_{\mathrm{min}}(z)=4\pi (H_0cz)^2f_{\mathrm{lim}}`$.

If the distribution of galaxy luminosities is described by the galaxy luminosity function $`\mathrm{\Phi }(L)`$, it is then simple to calculate what fraction of all galaxies is observable at any redshift:

$$\varphi (z)=\frac{_{L_{\mathrm{min}}(z)}^{\mathrm{}}𝑑L\mathrm{\Phi }(L)}{_0^{\mathrm{}}𝑑L\mathrm{\Phi }(L)}.$$
(1)
The quantity $`\varphi (z)`$ is usually referred to in the literature as the “selection function” peebles80a . The most common methods used to determine the galaxy luminosity function from the survey itself are those of efstathiou88a  and sandage79a . These methods assume that the luminosity function has universal shape, but not that galaxies are distributed homogeneously.

Figure 1 shows in the top panel the distribution of galaxies and redshifts in CfA2. Here we express galaxy luminosity in terms of the “absolute magnitude” $`M=2.5\mathrm{log}_{10}L+`$ const. The thick solid line shows the flux limit of the survey, translated into an absolute magnitude limit at each redshift. Because of this limit a number of galaxies which are observable at low redshifts are too faint to be observed at higher redshift.

The fraction of galaxies $`\varphi (z)`$ between absolute magnitudes $`22<M<10`$ which are unobservable at each redshift is shown in the bottom panel of Figure 1, based on a fit to the luminosity function in the survey using the method of efstathiou88a . Because the function falls rapidly from unity, it is clear that even at low redshifts the effects of the flux limit are important. We can use $`\varphi (z)`$ to calculate the expected distribution of observed galaxies with redshift. The top panel of Figure 2 compares these expected counts (dotted line) in redshift shells of thickness 0.001 to the observed counts (solid line) in CfA2. It appears that locally we are in an overdensity of galaxies of about a factor of two; note that at large distances, where each shell corresponds to a considerable amount of volume, the number of galaxies is very nearly the expected number.

If the flux limits are ignored, as in tanco , the incorrect conclusion that we live in a large overdensity is easily recovered (dashed lines in Fig. 2).

As mentioned above, the CfA2 survey covers a relatively small fraction of the sky. Thus, the PSC$`z`$ redshift survey provides a more useful sample to use in the context of this paper. Using the selection function provided by saunders00a , we again show the expected versus the observed counts for the PSC$`z`$ survey in the bottom panel of Figure 2. This survey also shows we are living in a slight overdensity, and furthermore reveals the general homogeneity of the nearby universe. (The actual counts and their dependence on redshift are slightly different than for CfA2, because the galaxies are selected in different ways).

## Calculation of the diffuse flux of UHECRs

We calculate the diffuse flux of UHECR protons numerically and compare our results with the analytical calculations carried out as in bg88 . Our propagation code includes pair production and photopion production as energy losses and also adiabatic energy losses due to the expansion of the universe. Since the inelasticity for pair production is very low, we consider it as a continuous energy loss process. The magnetic field is not included. A more detailed description of the numerical approach and of the results is reported in bbo00 .

In order to compare the results of the simulation with the observed statistics of events of operating or planned detectors, we generate the events in such a way that the total number of events above $`10^{19}`$ eV equals the observed number for the specific experiment under consideration.

The results of the code are checked versus the analytical results for the modification factor from single sources and from a diffuse distribution of sources. The agreement is excellent, and the effect of the fluctuations at energies larger than $`(34)\times 10^{19}`$ eV is evident (see below). On average the simulated flux is slightly larger than the analytical one, as expected on the basis of the stochasticity of the process of photopion production (on small distances there is an appreciable chance that some protons do not interact at all).

We investigate the effects of the real distribution of sources on the observed spectra of UHECRs, for different choices of the injection spectrum.

The possibility that a local overdensity of sources of UHECRs may help in solving the problem of the existence of events above the GZK cutoff goes back to beregri79  and is summarized in bible . The effect can easily be understood, since the severe photopion energy losses limit the maximum distance of UHECRs to distances of a few Mpc, while lower energy CRs can come from much larger distances. A local overdensity mainly affects, as a consequence, the fluxes of UHECRs. This is illustrated in Fig. 3 for a toy model in which the local overdensity is a top-hat function with $`\mathrm{\Delta }\rho /\rho =1,10,30`$ (solid, dashed and dash-dotted lines respectively) in a region of $`20`$ Mpc around the Earth. The fluxes of UHECRs have been calculated following the analytical approach of bg88 .

In order to have a direct comparison with the results of Ref. tanco , we first consider the case of an injection spectrum $`E^\gamma `$ with $`\gamma =3`$ and a source distribution extended up to a maximum redshift $`z_{max}=0.1`$, that corresponds approximately to the maximum redshift available in the PSCz catalog. It is worth stressing that the reason for the choice $`\gamma =3`$ in tanco  was motivated by the need to first reproduce the results of Yoshida and Teshima yoshida , whose curve is usually superimposed to the AGASA results AGASA . However, that curve seems to be obtained for $`\gamma =2.3`$ (Teshima, personal communication). The results of our calculations and of Ref. tanco  seem to be incompatible with such a choice of the power index.

We carry out the simulation for a homogeneous distribution of the sources and then for a source distribution that follows the profile found in the previous section, after the correction for selection effects.

The generation of events is ended when the total number of events with energy above some threshold equals the observed number. In fig. 4 we normalize the flux at $`10^{19}`$ eV and we stop the generation of new events when the total number of events above $`10^{19}`$ eV becomes 728 111 This sample includes the newly released AGASA data AGASAnew , equal to the number of events detected by AGASA in that energy range. The left upper panel is for a homogeneous distribution of the sources, while the right upper panel is obtained by adopting the PSCz distribution of galaxies. The diamonds are the AGASA data points and the crosses are the results of the simulation. The error bars in the simulation are obtained by generating 100 realizations and calculating their mean and variance. The continuous curves represent the result of the analytical calculation for the same value of the parameters. The lower curve is for the homogeneous case and the upper curve for the distribution derived from the PSCz catalog. One can easily see that the difference is small, which is expected since the correction for selection effects considerably reduces the local overdensity in comparison to what found by tanco . The total number of events with $`E>10^{20}`$ eV is $`1.2\pm 1.0`$ for the homogeneous case and $`1.5\pm 1.0`$ for the PSCz catalog galaxies. No one of our realizations had a number of events above $`10^{20}`$ eV comparable with observations. The main reason for that is that the spectrum is quite steep, in addition to the suppression due to photopion production.

In order to study the effect of an increased statistics of events, we simulated the situation for the case of Auger (9075 events predicted above $`10^{19}`$ eV in the first 3 years of operation). The situation is illustrated in the two lower panels of fig. 4 (on the left the homogeneous case and on the right the PSCz case). In general, the size of the error bars decreases everywhere. The number of events at $`E>10^{20}`$ eV (see figures) in the two cases is still much smaller than the projected one ($`100`$).

Since the main reason for the small number of events at high energy in the simulations is the steep spectra adopted above and in tanco , it is natural to look for different choices.

Many acceleration mechanisms produce spectra which are quite flatter than $`E^3`$. We calculate the expected flux of UHECRs for an injection spectrum $`E^{2.1}`$. In this case the normalization is chosen in such a way that the total number of events with energy $`E>4\times 10^{19}`$ eV is the same as observed by AGASA (49 events). The reason for this normalization will be clear below.

The results of the simulation and the corresponding analytical result are plotted in fig. 5 for three cases: homogeneous distribution of the sources with $`z_{max}=0.1`$ (crosses), PSCz distribution with $`z_{max}=0.1`$ (stars), homogeneous distribution with $`z_{max}=1`$ (squares).

The number of events at $`E>10^{20}`$ eV is mainly determined by the local distribution of the sources. For the adopted normalization, the homogeneous distribution gives $`3.3\pm 1.6`$ events above $`10^{20}`$ eV (to be compared with 8) and the PSCz distribution provides $`3.7\pm 2.0`$ events in the same range. In this last case, about $`5\%`$ of our realizations give a number of events above $`10^{20}`$ eV which is equal to or larger than the observed one. 222If the old AGASA statistics of 47 events above $`4\times 10^{19}`$ eV is used, $`15\%`$ of the realizations give a number of events above $`10^{20}`$ eV equal to or larger than the 6 observed then. The conclusion that the number of events above $`10^{20}`$ eV is within $`2`$ sigmas of the observations for an injection spectrum $`E^2`$ is consistent with the findings in bachall .

The deficit of events at energies lower than $`(34)\times 10^{19}`$ eV is evident in the case $`z_{max}=0.1`$. The high redshift sources only contribute additional flux at the low energies, therefore we also considered the case $`z_{max}=1`$. From fig. 5 it appears that the deficit is now evident only at energies lower than $`2\times 10^{19}`$ eV, where additional factors might further improve the low energy agreement. Some possible factors are: 1) cosmological magnetic fields; 2) source evolution; 3) a separate (possibly galactic) component which is relevant at lower energies and has a steep spectrum. These possibilities are discussed at length in bbo00 .

## Super Heavy Dark Matter: the direct connection

In this section we consider another example of a possible connection between the problem of the dark matter and the explanation of the UHECRs.

Heavy particles ($`m_X10^{12}10^{14}`$ GeV) can be produced in the Early Universe in different ways kolb ; bere ; kuzmin ; kt  and their lifetime can be finite though very long compared to the present age of the universe. In these circumstances, super-heavy particles can represent an appreciable fraction, if not all of the cold dark matter in the universe kolb ; bere ; kuzmin ; kt . The occasional decay of these particles results in the production of UHECRs, as widely discussed in the literature bere ; sarkar ; bbv ; blasi . In particular, if the relics cluster in galactic halos, as is expected, they can explain the cosmic ray observations above $`5\times 10^{19}`$ eV.

The decay of heavy relics results usually in the production of a quark-antiquark pair which rapidly hadronizes, generating two jets with approximately $`95\%`$ of the energy in pions, and $`5\%`$ in baryons. The decay of the pions results in the observed high energy particles, mainly in the form of gamma rays, and in the generation of ultra-high energy neutrinos. The spectrum of the gamma photons is relatively flat ($`E^{1.5}`$) reflecting the behaviour of the fragmentation function for the quarks. Therefore two main signatures of this model are: i) a flat energy spectrum; ii) composition dominated by gamma rays rather than by protons. Moreover, as in all top-down models, heavy elements are expected to be completely absent. Several calculations of the expected fluxes have been performed and presented in the literature. In Fig. 6 we report the results of the calculations of Ref. blasi .

The solid curves in fig. 6 are obtained by using a supersymmetric generalization of the quark fragmentation function, derived in SUSY , while the dashed lines are obtained by using the standard MLLA-QCD fragmentation function dok . The thick lines are for $`m_X=10^{14}`$ GeV, while the thin lines are for $`m_X=10^{13}`$ GeV.

In dt ; bbv  it was first recognized that, due to the asymmetric position of the Earth in the Galaxy, an appreciable anisotropy would result in this model. In bm ; medina  this issue was considered more quantitatively, taking into account the exposure of the present experiments. All authors concur that the present data is consistent with the predictions of the relic model for practically all reasonable values of the model parameters.

Recently, an interesting pattern has arisen from the analysis of the events with energy larger than $`4\times 10^{19}`$ eV: in AGASA  the sample with this energy cut comprises 47 events 333The data in AGASAnew  were not included in this analysis., whose overall distribution in space does not show appreciable deviation from isotropy. However, 3 doublets and one triplet were identified within an angular scale of $`2.5^o`$, comparable with the angular resolution of the experiment. A complete analysis, including the whole set of UHECR events above $`4\times 10^{19}`$ eV from the existing experiments was performed in watson . This extended sample comprises 92 events and shows 12 doublets and two triplets (each triplet is also counted as three doublets) within an angle of $`3^o`$. The chance probability of having more than this number of doublets was estimated to be $`1.5\%`$. Although it is probably too soon to rule out the possibility that these multiplets are just a random fluctuation, it is instructive to think about the possibility that their presence contains some physical information about the sources of UHECRs. Most of the top-down models for UHECRs (e.g. strings, necklaces, vortons, etc.) cannot naturally explain the multiplets.

Here we discuss how the multiplets can be interpreted in the context of the super-heavy dark matter (SHDM) model, following the discussion in bs00

We know a few things on dark matter, mainly as suggested by N-body simulations (see for instance simul ). Dark matter seems to be clustered in galactic halos with a distribution strongly peaked in the center. We model this distribution as in nfw :

$$n_H(r)=n^0\frac{(r/r_c)^1}{\left[1+\frac{r}{r_c}\right]^2},$$
(2)
where $`r_c`$ is the core size and $`n^0`$ is a normalization parameter. These two parameters can be set by requiring that the halo contains a given total mass ($`M_H`$) and that the velocity dispersion at some distance from the center is known (in the case of the Galaxy, the velocity dispersion is $`200`$ km/s in the vicinity of our solar system.). Alternative fits to the simulated dark matter halos and a discussion of whether or not simulated halos appear to be consistent with observations are provided in simul .

In addition to the smooth dark matter distribution, represented by eq. (2), N-body simulations also show that there is a clumped component which contains $`1020\%`$ of the total mass. The presence of these clumps are a natural consequence of the way in which gravity assembles dense virialized halos such as our galaxy today from the initially smooth density fluctuation field which was present when the cosmic microwave background (CMB) decoupled from the baryons. A more extended discussion on the formation and merging of the clumps can be found in bs00  and references therein.

We found that a good fit to the joint distribution in clump mass and position in the simulations of simul  is given by

$$n_{cl}(r,m)=n_{cl}^0\left(\frac{m}{M_H}\right)^\alpha \left[1+\left(\frac{r}{r_c^{cl}}\right)^2\right]^{3/2},$$
(3)
where $`n_{cl}^0`$ is a normalization constant, $`r_c^{cl}`$ is the core of the clumps distribution, and $`\alpha `$ describes the relative numbers of massive to less massive clumps. The simulations suggest that $`\alpha 1.9`$ simul . The constraints on the core size are weaker—we will study the range where $`r_c^{cl}`$ is between 3 and 30 percent of $`R_H`$. In simul , a halo with $`M_H2\times 10^{12}M_{}`$ contains about $`500`$ clumps with mass larger than $`10^8M_{}`$. This sets the normalization constant in eq. (3).

Clumps in the parent NFW halo are truncated at their tidal radii. The tidal radius of a clump depends on the clump mass, the density profile within the clump, and on how closely to the halo center it may have been. We assume that clumps of all mass are isothermal spheres (even though they are not truly isothermal, simul  suggest this is reasonably accurate): $`\rho _{cl}(r_{cl})1/r_{cl}^2`$, where $`r_{cl}`$ is the radial coordinate measured from the center of the clump.

As shown in bbv ; medina , the total (energy integrated) flux of UHECRs per unit solid angle from a smooth distribution of dark matter particles in the halo is $`\frac{d\mathrm{\Phi }}{d\mathrm{\Omega }}_0^{R_{max}}𝑑Rn_H(r(R))`$, where $`R`$ is the distance from the detector, and $`r`$ is the distance from the galactic center (so $`R`$ and $`r`$ are related by trigonometrical relations accounting for the off-center position of the Earth in the Galaxy). The upper limit, $`R_{max}`$, depends on the line of sight.

It is intuitively obvious that clumped regions will give an excess of events from certain directions, as was first pointed out in bereproc . Details of the numerical calculation of the fluxes and direction of arrivals of UHECRs can be found in bs00 .

Fig. 7 shows as an example one of the generated flux maps: the map represents the ratio of the total flux including the contribution from clumps, to the flux obtained by using a smooth NFW profile. The various free parameters were $`r_c=8`$ kpc, $`r_c^{cl}=10`$ kpc, and the mass distribution was truncated at a clump mass of $`1\%`$ of the mass of the NFW halo. This sort of plot emphasizes the clump contribution.

To calculate the small scale anisotropies, we generated $`10^4`$ mock samples, each of 92 observed events, and counted the number of doublets and triplets for angular scales of 3, 4, and 5 degrees. Our codes can also be used to check the corresponding numbers for the case of isotropic arrival directions (as in watson ). Two sets of values of the cores for the NFW and the clumped component were adopted, one in which $`r_c=8`$ kpc and $`r_c^{cl}=10`$ kpc (case 1) and the other with $`r_c=r_c^{cl}=20`$ kpc (case 2). The observed numbers of doublets within 3, 4, and 5 degrees for an isotropic distribution of arrival directions are given in watson  and are 12, 14 and 20 respectively. The number of doublets that we obtain in case 1 are 8, 14, and 21 within 3, 4, and 5 degrees respectively. The probability that the number of doublets equals or exceeds that observed is $`12\%`$, $`47\%`$ and $`57\%`$ respectively. This should be compared with the $`1.5\%`$, $`13.4\%`$ and $`15.9\%`$ quoted in watson  for an isotropic distribution of arrival directions.

We repeated the same calculation for the case 2. The corresponding averages and probabilities of exceeding the observed number of doublets within 3, 4 and 5 degree scales are 6.6, 12, and 18, and $`4.5\%`$, $`29\%`$ and $`36\%`$ respectively.

In both cases 1 and 2, the number of doublets on angular scales of 4 and 5 degrees is consistent with the observed values; presumably the discrepancy at 3 degrees is random chance.

We have also studied the occurrence of triplets. There is some ambiguity as to how a triplet is best defined; we have chosen to define triplets as configurations in which all three pairs would have been classified as doublets. (This means, for example, that a co-linear configuration of two doublets is not necessarily a triplet.) With this definition, the average number of triplets in case one is 0.5, 1.5 and 3, with the probability of having more than the observed triplets (2, 2, 3 respectively) equal to 4%, 16% and 35%. For case 2, the correspondent numbers are 0.4, 1, and 2.5 triplets and $`2\%`$, $`8\%`$ and $`20\%`$ for the probabilities to have more triplets than observed.

It is instructive to explore the reasons for the multiple events in the SHDM model. If we study the case in which all the halo mass is in the smooth NFW component, then the number of doublets typically drops by one or two. This suggests that the anisotropy due to our position in an NFW halo can result in a number of multiplets of events which is considerably larger than if the arrivals were from an isotropic background. The number of multiplets from the clumped component is mainly affected by the presence of large nearby clumps, whose number depends on the high mass cutoff imposed in the mass function of clumps. A maximum mass of $`1\%`$ of the halo mass implies a total mass in the clumps of $`1015\%`$ of $`M_H`$, consistent with the results of the simulations simul . Larger cutoffs imply larger mass fractions, which are harder to reconcile with the N-body simulations.

## Conclusions

We discussed here two possible connections between the problem of origin and clustering properties of dark matter and the nature of the particles with energy in excess of $`10^{20}`$ eV.

In the first part we investigated an indirect connection, consistent in the effect on the fluxes of UHECRs due to local clustering of the sources, driven by the corresponding clustering of dark matter. It has been known for a long time beregri79 ; bible  that a local overdensity in the sources of UHECRs may shift the energy of the GZK cutoff towards larger energies. The overdensities needed to make this effect relevant as a possible solution of the UHECR puzzle are of the order of $`10`$ or more. In a recent paper, Medina-Tanco tanco  claimed that this is exactly the value extracted from large scale structure surveys like the CfA2 catalog. We reconsidered this problem and showed that, by correctly extracting the density field from the CfA2 and the PSCz catalogs, the local overdensity is not larger than $`2`$, so that the calculations of the propagation of UHECRs are not particularly affected by the local source distribution. In other words the increase in the number of events above $`10^{20}`$ eV, in comparison with the case of homogeneous distribution is not important bbo00 .

On the other hand we noticed an interesting point: using an injection spectrum $`E^{2.1}`$ typical of astrophysical sources, we obtain a number of events with energy above $`10^{20}`$ eV, which is basically compatible with the observed number, once the statistical fluctuations and the Poisson noise in the photopion production have been taken into account. The agreement of the prediction at lower energies, depend on additional physics (evolution of the sources, cosmological magnetic fields or an additional (possibly) galactic component at energies around $`(12)\times 10^{19}`$ eV) that was not included in our calculations.

In the second part of the paper, we investigated a direct connection between dark matter and UHECRs. It has been proposed that most of the dark matter might be made of super heavy quasi-stable particles, created in the early universe. These particles would cluster in large scale structures, and in particular in the galactic halo. The rare decays of the super heavy dark matter particles in the halo can easily explain the fluxes of UHECRs that are observed. We reviewed the testable predictions of this model and then concentrated on a particular aspect, the anisotropy of arrival directions. A large scale (dipole) anisotropy in this model is easily foreseeable, since the earth is off center in the Galaxy. The predicted large scale anisotropy is compatible with the observed one bm ; medina , mainly due to a lack of exposure in the direction of the galactic center. Based on the clustering properties of dark matter on smaller scales, we also discussed the recent results of bs00  on the the small scale anisotropies. Based on the informations provided by N-body simulations of structure formation, the combination of the smooth peculiar density profile in the Galaxy and the clumped component survived in the halo can explain the observed small scale anisotropies, in terms of doublets and triplets of events. Future full sky experiments as Auger auger  will be extremely important to test this model, mainly in two ways: 1) an increased statistics of events will allow a better determination of the small angle clustering, and 2) a better composition determination will definitely allow to understand if gamma rays are an important component of the UHECRs.

Aknowledgments I am grateful to my collaborators, M. Blanton, A.V. Olinto and R.K. Sheth, for continuos interactions. I am also grateful to the organizers of the “International Workshop on Observing UHECRs from Space and Earth” for the nice environment that allowed instructive discussions with P. Biermann, J. Ellis, T. Gaisser, D. Harari, A. Letessier-Selvon, J. Linsley, L. Scarsi, G. Sigl, M. Teshima and T. Weyler among others. This work was supported by the DOE and the NASA grant NAG 5-7092 at Fermilab.


text:
# References

Energy calibration of large underwater detectors using stopping muons

V. A. Kudryavtsev, R. A. Brook, S. L. Cartwright, J. E. McMillan,
N. J. C. Spooner, L. F. Thompson

Department of Physics and Astronomy, University of Sheffield, Sheffield S3 7RH, UK

Abstract

We propose to use stopping cosmic-ray muons in the energy calibration of planned and deployed large underwater detectors. The method is based on the proportionality between the incident muon energy and the length of the muon path before it stops. Simultaneous measurements of the muon path and the amplitude of the signal from the photomultiplier tubes allow a relation between the energy deposited in the sensitive volume of the detector and the observed signal to be derived, and also provide a test of detector simulations. We describe the proposed method and present the results of simulations.

PACS: 95.55.Vj; 96.40.Tv; 14.60.Ef; 13.10.+q

Keywords: Muons underground; Stopping muons; Energy calibration; Underwater detectors; Neutrino telescopes; Monte Carlo simulations

Corresponding author: V. A. Kudryavtsev, Department of Physics and Astronomy, University of Sheffield, Hicks Building, Hounsfield Rd., Sheffield S3 7RH, UK

Tel: +44 (0)114 2224531; Fax: +44 (0)114 2728079;

E-mail: v.kudryavtsev@sheffield.ac.uk

1. Introduction

Several large neutrino telescopes have been proposed and are being deployed deep under water or under ice . They are designed to measure the fluxes and spectra of high-energy astrophysical neutrinos by detecting muons produced via neutrino-nucleus interactions. While these telescopes are designed to detect up-going muons produced by neutrinos passing right through the Earth, they can also detect down-going atmospheric muons. For most applications, such muons are a background which can be suppressed by several orders of magnitude by the measurements of the arrival times of Cherenkov photons and proper reconstruction of the muon track. These muons, however, being of known origin with predictable flux and spectrum, can serve also to test detector Monte Carlo and calibration techniques. One possible application of atmospheric muons is the energy calibration of the detector.

The energy of a muon above 1 TeV can be derived from the number of collected Cherenkov photons (converted to the amplitude or ADC channel). Since for high-energy muons the mean energy loss is proportional to the muon path and energy (stochastic energy losses dominate over ionisation loss), the muon energy loss inside the detector volume (number of produced Cherenkov photons) can be a measure of muon energy if the muon path is reconstructed. This concept of muon energy measurement is proven by various simulations . The simulations are usually based on full three-dimensional Monte Carlo simulations of the detector including muon energy loss, production and propagation of Cherenkov photons, response of photomultiplier tubes and electronics. The question is: how well can all the characteristics of all the elements of the detector be known or how accurately can they be measured with test light sources such as LEDs and lasers? Is it possible to check the simulation based predictions independently? In other words: can we know the muon energy from some other independent source, which we can then compare with Monte Carlo predictions and the measurement of the signal amplitude?

In the physics of keV-MeV energies there is a simple answer to a similar question. There are specific radioactive sources with mono-energetic X-ray or gamma-ray lines which are routinely used for energy calibration of detectors. Only accelerators can provide a quasi-monoenergetic particle beam at GeV-TeV energies, but they are usually located far from large underwater or underground detectors. Underground neutrino telescopes and cosmic-ray muon detectors, being smaller in size (compared to underwater detectors) and segmented structure, cannot measure the energy of high-energy muons but can use their mean or peak energy losses together with their path to calibrate the detectors (see, for example, ). For the larger underwater detectors this method can hardly be applied because the peak in the muon energy loss distribution is smoothed by stochastic losses due to bremsstrahlung, pair production and inelastic scattering over the long muon path.

We propose to use stopping atmospheric muons for this task. They have the advantage that their incident energy can, in principle, be measured through their path in the sensitive volume of the detector and, hence, can be compared to the observed signal and predictions of simulations. In Section 2 the basic idea of calibration is described and in Section 3 main requirements for the detector hardware and software are given. Results of our simulations for a simplified detector are presented in Section 4 and the conclusions are given in Section 5.

2. Energy calibration of underwater detectors

Usually energy calibration means the relation between particle energy and the amplitude of the signal processed by DAQ. This definition is correct if the particle energy is totally absorbed in the detector. For underground detectors most muons are not absorbed in the detector and the relation can be derived only between the observed signal and the energy deposited in the detector. Usually it is assumed that this relation is linear.

Since there is no mono-energetic source of muons underwater, normal methods of calibration cannot be applied. A mono-energetic source is not necessary, however, if we can determine the muon energy from another measurement. There is no way to measure the energy of individual muons from the whole muon population, but there exists a sub-sample consisting of stopping muons for which the incident energy can be known from the measurement of the muon path length between the point where the muon enters the sensitive volume of the detector and the point where it stops. The idea of the method is to measure the muon path if it stops within the sensitive volume of the detector, to calculate its initial energy using the energy-range relation and to compare it to the measured signal amplitude, thus deriving the relation between the amplitude and the deposited energy and checking Monte Carlo predictions. Subsequently, the relation between the measured amplitude and the energy deposition can be used to estimate the energy of through-going TeV muons.

One of the crucial points of the method is an estimation of incident muon energy from the muon path length. Simple considerations show that, since low energy muons (only such muons stop in the detector) mainly lose energy via ionisation and atomic excitation without undergoing stochastic energy loss, for most muons a simple linear relation will work without large fluctuations. The results of the simulations will be shown in Section 4.

Another important question is: can large underwater detectors measure path length of a stopping muon with sufficient accuracy? This is discussed in the following section.

3. Detector performance

What are the performance requirements of an underwater detector for it to be able to use stopping muons for energy calibration?

1. Such a detector should be sensitive to down-going atmospheric muons. This is not necessarily true since neutrino telescopes have most (and possibly all) optical modules looking downwards. Measurements of the depth-intensity relation for down-going muons by the AMANDA and Baikal Collaborations  prove that down-going muons can be reconstructed and their fluxes can be measured.

2. Separation between optical modules should not be large to allow i) the detection of weak signals from low-energy muons and ii) the determination of the stopping point of the muon with sufficient accuracy. Requirement i) may reduce the sensitive volume of the detector for this task but is unlikely to make the energy calibration impossible. (Note, that the Baikal energy threshold has been estimated as 10 GeV  which is well matched to our requirements). The requirement ii) is the most important. The distance between optical modules on one string varies from 6 to 20 metres, while the distance between strings is 20-80 metres for different arrays. Simulations performed by the ANTARES Collaboration  for neutrino oscillation studies show, however, that contained and semi-contained neutrino-induced events can be successfully discriminated from through-going muons and, hence, both the points of muon production and absorption can be determined. Full reconstruction of muon events is more difficult, of course, in the case of down-going muons because of their large number. But even a fraction of stopping muon events (well-reconstructed) may be enough for energy calibration. (Note, however, that there is a danger of selecting a sub-sample of reconstructed stopping muons with characteristics different from those of the whole population of stopping muons).

3. The rate of stopping muon events should be large enough to perform calibration in a feasible period of time. This requirement may be more difficult to satisfy at deep sites. We will present our estimates in Section 4.

4. Even if all requirements are satisfied, the proposed method cannot avoid the use of detector simulation. This is due to the dependence of the observed amplitude on the distance to the muon track. Imagine two vertical muon tracks, one close to the string of optical modules and another far away from the string. Even if the lengths of the tracks are equal and the incident muon energies are equal too, the observed signals will be different due to different distances from the tracks to the phototubes. Such simulations have already been performed by the collaborations . The results depend strongly on the properties of water/ice. The propagation of Cherenkov photons is not included in our simulations, the assumption being that appropriate corrections for the distance between optical modules and muon track can be made for a particular detector if the track coordinates are fully reconstructed.

5. Finally, multiple muon events can present a source of background for stopping muons. Again, measured depth-intensity curves  prove that at least some of the muon tracks can be reconstructed properly.

As previously mentioned, there is a problem of possible bias in the calibration if the sample of stopping muons that can be reconstructed has particular characteristics different from those of the whole population (for example, energy loss higher than the mean value). However, the probability of enhanced energy loss along the whole muon track is extremely low. High amplitude photomultiplier hits due to occasional muon-induced local cascades can be excluded from the analysis. The presence of bias can be detected by comparing the shape of the measured distribution of specific energy loss with the simulated distribution.

4. Simulation of stopping muons

In this section the results of simulations are presented which demonstrate that energy calibration, in principle, can be performed using stopping muons. Further investigation of this question can be made for each individual detector using appropriate software packages including the detector design and characteristics.

The simulation was done in the following way: muons were sampled according to their energy spectrum and angular distribution at a particular depth underwater, propagated through the detector volume and their initial and final characteristics were stored on disk. Energy spectra and angular distribution of muons under water were obtained by propagating muons with various initial energies from sea level down to various depths in water using the muon transport package MUSIC . Then, muon energy distributions underwater were convolved with muon spectra at sea level using the relevant parameterisation from .

To check the correctness of the procedure we calculated the muon “depth – vertical intensity” relation and compared it with those measured by the Baikal  and AMANDA  experiments. The results are shown in Figure 1. Measured relations are well fitted to the calculated relation assuming the power index of the primary spectrum is equal to 2.78 – in good agreement with underground experiments .

A special package was developed to sample single muons according to their spectra and angular distribution underwater: MUSUN (MUon Simulation UNderwater/UNderground) . As an example, we present here the results for the depth of 2 km under water. This corresponds approximately to the depth at which ANTARES will be deployed  and is 0.5 km deeper than the AMANDA site. A simple detector in the form of cylinder with a radius of 100 metres and a height of 500 metres was used. This cylinder represented the sensitive volume of the detector. Muons were sampled uniformly on the upper surface of the cylinder. We required also that the muon track (real or extended if the muon stopped) crossed the lower flat surface of the cylinder. It was assumed that the muon track and stopping point can be accurately reconstructed if the muon traversed between 100 metres and 400 metres of water inside the cylinder.

Figure 2 shows a scatter plot of initial energy versus track length for stopping muons with track lengths between 100 and 400 metres. No systematic uncertainty in track reconstruction was included. There is a clear proportionality between track length and initial energy with small fluctuations for most muons. Some of them, however, undergo, a large stochastic energy loss. Figure 2 can be converted into the distribution of specific energy losses (per unit track length). This is shown by a dash-dotted line in Figure 3 for the same stopping muon sample. The distribution is characterised by a narrow prominent peak which corresponds to the most probable specific energy loss. In practice, the track length can be reconstructed with some finite accuracy. If we assume that the reconstructed track length distribution follows a Gaussian distribution with mean value equal to the true value of the track length and standard deviation equal to 40 metres (at least twice the vertical distance between the modules on the string) the peak becomes smoother but is still present (dashed line in Figure 3). (Note, that the simulations of the ANTARES performance for neutrino oscillation study  show that the starting and stopping points of more than half up-going muons can be reconstructed with an accuracy of better than 30 m).

All previous results have been obtained assuming no fluctuations in the number of detected photons due to various distances between the optical modules and muon track. Such assumption is equivalent to the precise reconstruction of the position of the muon track with respect to the optical modules and correction for such fluctuations. It is unlikely, however, that this scenario is realised in practice. To account for the fluctuations in the number of collected photons in our simulations we sampled the reconstructed energy deposition according to a Gaussian distribution on the log(E) scale with mean value equal to the true energy deposition (which came as an output from the muon propagation code) and standard deviation equal to 0.176 on the log(E) scale (factor 1.5). The choice of the standard deviation for simulations was based on the known fact that the energy of multi-TeV through-going muons can be reconstructed with an accuracy of 0.4-0.5 on the log(E) scale . For high-energy through-going muons, however, the energy resolution is determined mainly by the energy loss fluctuations, while this effect (being much smaller for low-energy stopping muons) is included separately in our simulations. We included also the statistical fluctuations in the number of detected photoelectrons assuming that on average an atmospheric muon produces about 1 photoelectron per 1 GeV of energy deposition. To derive this number we used the mean number of photoelectrons produced by an atmospheric muon in the Baikal detector, quoted in  (23 photoelectrons), the length of Baikal string which is equal to 72 metres, and the mean energy deposition in water which is about 300 MeV/m. We assumed also that the energy threshold of the detector is equal to 20 GeV. The results of the simulations are shown in Figure 3 by solid curve. The distribution of specific energy losses is broader than in the simpler case when we did not account for the fluctuations in the number of detected photons. For comparison, we present also a similar distribution for through-going muons which pass more than 400 metres of water in the sensitive volume of the detector. It is obvious, that the distribution of specific energy losses for stopping muons is more favourable for the use in the energy calibration procedure. This does not exclude, however, the possibility to use such a distribution for the through-going muon sample either as an additional test or in the case that our assumptions about possible fluctuations were too optimistic.

The position of the peak calculated for each individual detector can be used to convert the measured track length of a stopping muon to its initial energy which can then be compared to the measured signal. For our simple detector at 2 km depth the rate of stopping muons which can be reconstructed is about 0.3 s<sup>-1</sup>.

Another important issue which we want to address in the discussion, is the contamination from multiple muons. Multiple muons contribute from a few percent to a few ten percent of the total muon flux even for very large detectors. Their contribution decreases with depth. Such events produce more light and hence, may be detected with higher probability than single muons. If so, they can smooth the observed distribution of specific energy losses.

To estimate the effect from multiple muons we simulated the development of vertical Extensive Air Showers initiated by primary protons using the code CORSIKA . Then, we propagated all muons down to 2 km of water with MUSIC and repeated simulations of muons (both single and multiple) in the detector. The results are presented in Figure 4. Solid curve shows the distribution of specific energy losses for single muons. Note that the distributions shown by solid curves in Figures 3 and 4 are very similar, though only vertical muons have been simulated with CORSIKA. Then, we assumed that only events with total observed energy more than 20 GeV, regardless of the muon multiplicity, can be detected. We assumed also that multiple muon events could be reconstructed as stopping single muons only if all muons stopped inside the detector. In this case the reconstructed track length corresponded to the muon with the longest range. The results are shown in Figure 4 by dashed curve. From Figure 4 we conclude that multiple muons do not change significantly the distribution of specific energy losses for stopping muons. The contribution of multiple muons from heavy primaries and from inclined directions can modify slightly calculated distribution of specific energy losses but is is unlikely to change this conclusion.

It is possible that the number of Cherenkov photons emitted by a low-energy muon is not enough to hit a phototube if the muon does not pass very close to it. In this case, only muons passing very close to the strings can be used for calibration and a single string analysis technique can be used (each string represents a separate detector). The basic features of the method remain the same, however, due to the reduced acceptance, the rate of reconstructed stopping muons will be approximately 2 orders of magnitude less (if a 10 m radius single string detector is used).

5. Conclusions

A method of energy calibration of large underwater detectors using stopping muons has been described. The method is proven to work by Monte Carlo simulations for a simple (“ideal”) detector. Further investigations involving structure and characteristics of specific neutrino telescopes are needed.

6. Acknowledgements

The authors wish to thank PPARC for financial support. The idea of the method and first preliminary results of the simulations were reported in May of 1998 on one of the meetings of the ANTARES Collaboration. The authors are grateful to the Collaboration for useful comments and discussion. One of the authors (VAK) thanks Dr. E. V. Korolkova for discussion and remarks. We wish to thank also anonymous referee for useful comments.


text:
# Lifetimes of the 𝑏–flavored baryons in the light–front quark model The talk presented at the 4th International Conference on Hyperons, Charm and Beauty Hadrons, Valencia, June 27–30, 2000. This work was supported in part by the RFBR grants, ref. No 00–02–16363 and 00-15–96786

## 1 INTRODUCTION

The lifetimes of the $`b`$ flavored hadrons $`H_b`$ are related both to the CKM matrix elements $`|V_{cb}|`$ and $`|V_{ub}|`$ and to dynamics of $`H_b`$ decays. Weak decays of heavy hadrons involve three fundamental scales, the weak interaction scale $`M_W`$, the $`b`$–quark mass $`m_b`$, and the QCD scale $`\mathrm{\Lambda }_{QCD}`$, which are strongly ordered: $`M_Wm_b\mathrm{\Lambda }_{QCD}`$. The first inequality means that the decay width of the $`b`$ quark rises so rapidly with $`m_b`$, namely $`\mathrm{\Gamma }_bG_F^2m_b^5`$, that the spectator process dominates for large $`m_b`$. The second inequality means that for large $`m_b`$ $`b`$–quark decays before it can hadronize. Thus in the limit $`m_b\mathrm{}`$ the lifetimes of all $`b`$ hadrons must be equal: $`\tau _{B^+}=\tau _{B_d^0}=\tau _{B_s^0}=\tau _{\mathrm{\Lambda }_b}`$.

Inclusive $`H_b`$ decays can be treated with the help of an operator product expansion (OPE) combined with the heavy quark expansion . The OPE approach predicts that all corrections to the leading QCD improved parton terms appear at the order $`1/m_b^2`$ and beyond. Thus mesons and baryons containing b quark are expected to have lifetimes differing by no more than a few per cent. The result of this approach for the $`\mathrm{\Lambda }_b`$ lifetime is puzzling because it predicts that  $`\left(\tau _{\mathrm{\Lambda }_b}/\tau _B\right)_{OPE}=0.98+𝒪(1/m_b^3)`$, whereas the experimental findings suggest a very much reduced fraction

$$\left(\tau _{\mathrm{\Lambda }_b}/\tau _B\right)_{exp}=0.78\pm 0.04,$$
(1)
or conversely a very much enhanced decay rate.

Application of OPE to decays of b–hadrons has potentially two caveats. One is that OPE is used in Minkowsky domain, and therefore relies on assumption of quark–hadron duality at the energies involved in decays. Another uncertainty arises from poor knowledge of matrix elements of four–quark operators arising as terms of OPE. The decay rates of B and $`\mathrm{\Lambda }_b`$ are $`\mathrm{\Gamma }(B)=0.63\pm 0.02\mathrm{ps}^1`$ and $`\mathrm{\Gamma }(\mathrm{\Lambda }_b)=0.83\pm 0.05\mathrm{ps}^1`$ differing by $`\mathrm{\Delta }\mathrm{\Gamma }(\mathrm{\Lambda }_b)=0.20\pm 0.05\mathrm{ps}^1`$. The four–fermion processes of weak scattering and Pauli interference calculated using the factorization approach and the description of the baryon relying on quantum mechanics of only the constituent quarks could explain, under certain conditions, only $`(13\pm 7)\%`$ of $`\mathrm{\Delta }\mathrm{\Gamma }(\mathrm{\Lambda }_b)`$. . However, the ”spectator effects”, which formally appear at $`𝒪(1/m_b^3)`$ in OPE can be considerably larger. A piloting lattice study of Ref.  suggests that the effects of weak scattering and interference can eliminate $`50\%`$ of the discrepancy between the theoretical prediction for $`\tau _{\mathrm{\Lambda }_b}/\tau _B`$ and the experimental result. It also appears that not all of the discrepancy can be accounted for by spectator effects. The OPE prediction, given all the uncertainties involved, is unlikely to produce a number lower that $`0.9`$.

It is not clear whether the present contradiction between the theory and the data on $`\tau _{\mathrm{\Lambda }_b}/\tau _B`$ is a temporary difficulty or an evidence of fundamental flaws in theoretical understanding. In spite of great efforts of experimental activity the $`\mathrm{\Lambda }_b`$ lifetime remains significantly low which continues to spur theoretical activity. In this respect, the use of phenomenological models, like the constituent quark model, could be of interest as a complimentary approach to the OPE resummation method. In this talk I shall consider the preasymptotic effects for the $`\mathrm{\Lambda }_b`$ lifetime in the framework of the light–front (LF) quark model. My talk is preceded by that of Dr.Kalman , covering, in part, the same topic. Therefore, I will avoid those aspects that have been already discussed. Instead, I will concentrate on the calculation of lifetimes of $`\mathrm{\Lambda }_b`$ and other $`b`$–flavored baryons.

## 2 LF QUARK MODEL

The LF quark model  is a relativistic constituent quark model based on the LF formalism. In Ref.  this formalism has been used to establish a simple quantum mechanical relation between the inclusive decay rate of the B meson and that of a free $`\mathrm{b}`$ quark. The approach of  relies on the idea of duality in summing over the final hadronic states. It has been assumed that the sum over all possible charm final states $`X_c`$ can be modeled by the decay width of an on–shell $`b`$ quark into on–shell $`c`$ quark folded with the $`b`$–quark distribution function $`f_B^b(x,p_{}^2)=|\psi _B^b(x,p_{}^2)|^2`$. The latter represents the probability to find $`b`$ quark carrying a LF fraction $`x`$ of the hadron momentum and a transverse relative momentum squared $`p_{}^2`$. For the partial rates $`d\mathrm{\Gamma }(B)/dq^2`$ the above mentioned relation takes the form

$$\frac{d\mathrm{\Gamma }(B)}{dq^2}=\frac{d\mathrm{\Gamma }^b}{dq^2}R_B(q^2),$$
(2)
where $`d\mathrm{\Gamma }^b/dq^2`$ is the free quark semileptonic or nonleptonic differential decay rate, with $`q`$ being the 4–momentum of the $`W`$ boson. The function $`R_B(q^2)`$ incorporates the non perturbative effects related to binding effect and to the Fermi motion of the heavy quark inside the hadron and is expressed in terms of the bound–state factor $`|\psi _B^b(x,p_{}^2)|^2`$. For the details see . The non–perturbative corrections to the free quark decay rate vanish as $`m_b\mathrm{}`$  but may be essential at finite values of the $`b`$ quark mass due to the difference between $`m_b`$ and $`m_{H_b}`$ and primordial motion of the $`b`$ quark inside $`H_b`$.

To define the LF wave function we use the simple operational anzätz  and express $`\psi _B^b(x,p_{}^2)`$ in terms of the equal time radial wave function $`\psi _B(p^2)`$. Explicit formulae relating $`\psi _B^b(x,p_{}^2)`$ and $`\psi _B(p^2)`$ can be found e.g. in Ref. .

## 3 $`B`$–MESON LIFETIME

In what follows the B meson orbital wave function is assumed to be the Gaussian function as $`\psi _B(p^2)=\left(1/\beta _{b\overline{d}}\sqrt{\pi }\right)^{\frac{3}{2}}\mathrm{exp}\left(p^2/2\beta _{b\overline{d}}^2\right)`$, where the parameter $`1/\beta _{b\overline{d}}`$ defines the confinement scale. We take $`\beta _{b\overline{d}}=0.45\mathrm{GeV}`$. For $`|\mathrm{\Psi }_B(0)|^2`$, the square of the wave function at the origin, we have $`|\mathrm{\Psi }_B(0)|^2=1.64\times 10^2\mathrm{GeV}^3`$ that compares favorably with the estimation in the constituent quark ansätz  $`|\mathrm{\Psi }_B(0)|^2=M_Bf_B^2/12=(1.6\pm 0.7)\times 10^2\mathrm{GeV}^3`$ for $`f_B=190\mathrm{MeV}`$.

We have found that the $`\tau _{\mathrm{\Lambda }_b}/\tau _B`$ ratio is rather stable with respect to the precise values of the heavy quark masses $`m_b`$ and $`m_c`$ provided $`m_bm_c3.5\mathrm{GeV}`$. From now on we shall use the reference values $`m_b=5.1\mathrm{GeV}`$ and $`m_c=1.5\mathrm{GeV}`$. The value of the CKM parameter $`|V_{cb}|`$ cancels in the ratio $`\tau _{\mathrm{\Lambda }_b}/\tau _B`$, but is important for the absolute rates. Details of our calculations of $`\mathrm{\Gamma }(B)`$ are given in Table 1 for the three values of the constituent mass $`m_{sp}`$. The values $`m_{sp}300(200)\mathrm{MeV}`$ are usually used in non-relativistic (relativistic) quark models. We have also considered a very low constituent quark mass $`m_{sp}=100\mathrm{MeV}`$ to see how much we can push up the theoretical prediction for $`\mathrm{\Gamma }(\mathrm{\Lambda }_b)/\mathrm{\Gamma }(B)`$, see below. All the semileptonic widths include the pQCD correction as an overall reduction factor equal to 0.9. Following Ref.  the transitions to baryon-antibaryon $`(\mathrm{\Lambda }_c\overline{N}\mathrm{and}\mathrm{\Xi }_{cs}\overline{\mathrm{\Lambda }})`$ pairs are included. In addition we have added $`\mathrm{BR}1.5\%`$ for the $`bu`$ decays with $`|V_{ub}/V_{cb}|0.1`$. The value of $`|V_{cb}|`$ is defined by the condition that the calculated B lifetime is $`1.56ps`$.

## 4 $`\mathrm{\Lambda }_b`$ LIFETIME

We shall analyze the inclusive semileptonic and non–leptonic $`\mathrm{\Lambda }_b`$ rates on simplifying assumption that $`\mathrm{\Lambda }_b`$ is composed of a heavy quark and a light scalar diquark with the effective mass $`m_{ud}`$. Then the treatment of the inclusive $`\mathrm{\Lambda }_b`$ decays is simplified to a great extent and one can apply the model considered above with the minor modifications. For the heavy–light diquark wave function $`\psi _{\mathrm{\Lambda }_b}`$ we again assume the Gaussian ansätz with the oscillator parameter $`\beta _{bu}`$. The width of $`\mathrm{\Lambda }_b`$ can be obtained from that of $`B`$ by the replacements $`M_BM_{\mathrm{\Lambda }_b}`$, $`m_{sp}m_{ud}`$ and $`\beta _{b\overline{d}}\beta _{bu}`$. The latter two replacements change $`f_{\mathrm{\Lambda }_b}^b`$, the $`b`$ quark distribution function inside the $`\mathrm{\Lambda }_b`$, in comparison with $`f_B^b`$.

The inclusive nonleptonic channels for $`\mathrm{\Lambda }_b`$ are the same as for B meson except the decays into baryon–antibaryon pairs which are missing in case of $`\mathrm{\Lambda }_b`$. The absence of this decay channel leads to the reduction of $`\mathrm{\Gamma }(\mathrm{\Lambda }_b)`$ by $`7\%`$. This effects prolongs $`\tau _{\mathrm{\Lambda }_b}`$ over $`\tau _B`$. The phase space enhancement in $`\mathrm{\Lambda }_b`$ is marginal and can not be responsible for a shorter lifetime of $`\mathrm{\Lambda }_b`$. The only distinction between the two lifetimes, $`\tau _{\mathrm{\Lambda }_b}`$ and $`\tau _B`$, can occur due to the difference of the binding and Fermi motion effects. These effects are encoded in the parameter $`x_0=m_b/m_{H_b}`$ and the distributions $`f_{\mathrm{\Lambda }_b}^b`$ and $`f_B^b`$.

The concrete amount of the deviation between $`\tau _{\mathrm{\Lambda }_B}`$ and $`\tau _B`$ depends in essential way on $`m_{ud}`$ and $`\beta _{ub}`$. In what follows the diquark mass is taken as

$$m_{ud}=m_{}=\frac{1}{2}(m_u+m_dm_\pi ).$$
(3)
In this relation inspired by the quark model, the factor $`1/2`$ arises from the different color factors for $`u`$ and $`\overline{d}`$ in the $`\pi `$–meson ( a triplet and antitriplet making a singlet) and $`u`$ and $`d`$ in the the $`\mathrm{\Lambda }_b`$ (two triplets making an antitriplet).

The quantity $`\beta _{ub}`$ can be translated into the ratio of the squares of the wave functions determining the probability to find a light quark at the location of the $`b`$ quark inside the $`\mathrm{\Lambda }_b`$ baryon and $`B`$ meson, i.e. $`r=|\mathrm{\Psi }_{\mathrm{\Lambda }_b}(0)|^2/|\mathrm{\Psi }_B(0)|^2=\left(\beta _{bu}/\beta _{b\overline{d}}\right)^3`$. If we take $`r=1.2\pm 0.2`$ as suggested by the lattice study of Ref.  we obtain $`\left(\tau _{\mathrm{\Lambda }_b}/\tau _B\right)_{LF}=0.88\pm 0.02`$ for $`m_{sp}=200\mathrm{MeV}`$. Adding the contribution of the 4–quark operators $`\mathrm{\Delta }\mathrm{\Gamma }^{4q}=0.08\pm 0.02`$ from  we reproduce the experimental result (1).

It is instructive to note that the calculated branching fractions of $`\mathrm{\Lambda }_b`$ show marginal dependence on the choice of the model parameters; they are $`11.5\%`$ for the semileptonic $`bce\nu _e`$ transitions, $`2.8\%`$ for $`bc\tau \nu _\tau `$, $`50\%`$ for the nonleptonic $`bcd\overline{u}`$ transitions, and $`16\%`$ for $`bc\overline{c}s`$ transitions.

An exploratory study of the dependence of $`\tau _{\mathrm{\Lambda }_b}`$ on $`m_{ud}`$ and $`r`$ has been performed in Ref. . Two choices for diquark masses were taken: $`m_{ud}=m_u+m_d`$ corresponding to zero binding approximation and $`m_{ud}=m_{}`$ from (3). The wave function ratio has been varied in the range $`0.3r2.3`$ 111 Recall that estimates of the parameter $`r`$ using the non relativistic quark model or the bag model or QCD sum rules are typically in the range $`0.10.5`$. The value $`r=2.3`$ corresponds to the upper Rosner estimation  obtained from the ratio of hyperfine splitting between $`\mathrm{\Sigma }_b`$ and $`\mathrm{\Sigma }_b^{}`$ baryons and $`B`$ and $`B^{}`$ mesons . that corresponds to $`0.3\mathrm{GeV}\beta _{bu}0.6\mathrm{GeV}`$. The results for $`\tau _{\mathrm{\Lambda }_b}/\tau _B`$ are exhibited in Table 2, which we now discuss.

We observe that to decrease the theoretical prediction for $`\tau _{\mathrm{\Lambda }_b}`$ requires to decrease the value of the hadronic parameter $`r`$ to 0.3-0.5 and the value of $`m_{sp}`$ to $`100\mathrm{MeV}`$. For example, assuming that $`r0.3`$ we find that the lifetime ratio is decreased from $`0.88`$ to $`0.81`$ if $`m_{sp}`$ is reduced from $`300\mathrm{MeV}`$ to $`100\mathrm{MeV}`$ and $`m_{ud}=m_u+m_d`$. For the $`m_{ud}=m_{}`$ the ratio is almost stable ($`0.8`$), so that reducing of the diquark mass produces a decrease of the lifetime ratio by $`1\%`$, $`5\%`$, and $`8\%`$ for $`m_{sp}=100,200,\mathrm{and}300\mathrm{MeV}`$, respectively. Varying the spectator quark mass in a similar way we find that for the ”central value” $`r1`$ the lifetime ratios are reduced from $`0.93`$ to 0.88 for $`m_{ud}=m_u+m_d`$ and remain almost stable ($`0.86`$) for $`m_{ud}m_{}`$. For the largest possible value of $`r`$ suggested in , $`r2.3`$, the lifetime ratios are reduced from 0.97 to 0.94 in the former case and remain almost stable $`0.910.93`$ in the latter case.

The lifetimes of other $`b`$–flavored baryons are reported in Table 3. The results are obtained using $`r=1`$ and adding the marginal contribution of $`\mathrm{\Delta }\mathrm{\Gamma }^{4q}`$ from Ref. . As we have discussed earlier this contribution may be strongly enhanced by the QCD effects in which case the lifetime predictions should be decreased by $`10\%`$.

## 5 CONCLUSIONS

If the current value of $`\left(\tau (\mathrm{\Lambda }_b)\right)_{exp}`$ persists, the most likely explanation is that some hadronic matrix elements of four–quark operators are larger than the naive expectation . If a significant fraction $`50\%`$ of the discrepancy between the theoretical prediction for $`\tau _{\mathrm{\Lambda }_b}/\tau _B`$ and the experimental result can be accounted for the spectator effects  then the remainder of the discrepancy can be naturally explained by the preasymptotic effects. We have found that the binding effect produces for $`r=1.2\pm 0.2`$ $`m_{sp}=200\mathrm{MeV}`$ and $`m_{ud}=250\mathrm{MeV}`$ an additional increase of $`\mathrm{\Delta }\mathrm{\Gamma }(\mathrm{\Lambda }_b)/\mathrm{\Gamma }(B)`$ by $`12\pm 2\%`$ 222Note that for these values of quark and diquark masses the Fermi motion effects in $`\mathrm{\Lambda }_b`$ and $`B`$ do not differ significantly. .

## ACKNOWLEDGMENTS

I would like to thank Calvin Kalman, Pepe Salt and Miguel–Angel Sanchis–Lozano for organizing an excellent Conference with a stimulating scientific program.