text:
# The Cosmic-Ray antiproton flux between 3 and 49 GeV.

## 1 Introduction

Detailed measurements of the cosmic-ray antiproton energy spectrum provide important information concerning the origin and propagation of cosmic rays in the Galaxy. In fact, antiprotons are a natural product of interactions between cosmic rays and the interstellar matter. Moreover, antiprotons can be produced by exotic sources such as evaporation of primordial black holes (Hawking, 1974; Kiraly, Wdowczyk, & Wolfendale, 1981; Maki, Mitsui, & Orito, 1996) and annihilation of supersymmetric particles (Stecker, Rudaz, & Walsh, 1985; Bottino et al., 1998; Bergström, Edsjö, & Ullio, 1999a, b). The measurement of the antiproton spectrum at energies above a few GeV permits the study of these topics free of uncertainties associated with the secondary antiproton production such as nuclear sub-threshold effects and of the uncertainties in the solar modulation effect. Furthermore, it permits us to investigate the possible annihilation of heavy supersymmetric particles (Ullio, 1999).

Several measurements of the cosmic-ray antiprotons have been performed since the first detection by Golden et al. (1979). However, most of these experiments have measured the antiproton spectrum at energies below 4 GeV (see Orito et al. (2000) and references therein). Only two experiments (Golden et al., 1979; Hof et al., 1996) have obtained antiproton results at energies above 4 GeV and these results differ by a large amount. We report in this paper a new observation of antiprotons with energies up to 50 GeV obtained with the CAPRICE98 experiment. This apparatus was launched by balloon from Fort Sumner, New Mexico (34.3 North Latitude, 110.1 West Longitude) and landed close to Holbroke, Arizona (34. North Latitude, 104.1 West Longitude), on May 28 and 29, 1998, at an atmospheric pressure of 4.2 to 6.2 mbar for 21 hours and average value of the vertical cut-off rigidity of about 4.3 GV (Shea & Smart, 1983).

Preliminary results on the antiproton to proton ratio from CAPRICE98 were reported earlier (Bergström et al., 2000). Here, we present the absolute energy spectrum of antiprotons in the energy region at the top of the atmosphere between 3 and 49 GeV. We also describe in detail the analysis of the flight data. The detector system is described in Section 2, the data analysis in Section 3 and the results are presented and discussed in Section 4.

## 2 The CAPRICE98 apparatus

Figure 1 shows the NMSU-WiZard/CAPRICE98 spectrometer (Ambriola et al., 1999). It included from top to bottom: a gas Ring Imaging Cherenkov (RICH) detector, a time-of-flight (ToF) system, a tracking system consisting of drift chambers and a superconducting magnet and a silicon-tungsten imaging calorimeter.

### 2.1 The Gas RICH detector

The Ring Imaging Cherenkov (RICH) detector was designed primarily to identify antiprotons in the cosmic rays in a large background of electrons, muons and pions (Carlson et al., 1994). The RICH detector (Francke et al., 1999; Bergström et al., 2001) consisted of a photosensitive multiwire proportional chamber (MWPC) and a 1 m tall radiator box, filled with high purity C4F10 gas. When a charged particle with $`\beta >1/\mathrm{n}`$, where n is the refractive index and $`\beta =v/c`$ ($`v`$ being the particle velocity and $`c`$ the speed of light), passed through the RICH detector, Cherenkov light was emitted in the gas radiator along the trajectory. The Cherenkov light was emitted at an angle ($`\theta _c`$) determined by the Cherenkov relation ($`cos(\theta _c)=\frac{1}{\beta \mathrm{n}}`$), and creating a cone of light in the same direction as the particle trajectory.

The cone of light after traversing the radiator volume was reflected back and focused by a spherical mirror toward the MWPC. There the Cherenkov light interacted with a photosensitive gas, tetrakis-dimethyl-amino-ethylene (also called TMAE), and photoelectrons were produced. These electrons were amplified and then detected by induced pulses in a matrix pad plane. This plane had an area of 51.2$`\times `$51.2 cm2, divided in $`64\times 64`$ pads of size $`8\times 8`$ mm2, where the cone of Cherenkov light gave a ring-like image. The size of the ring was dependent on the velocity of the particle. The ring diameter increased from 0 at the RICH threshold (about 18 GV for protons) to about 11 cm for a $`\beta 1`$ particle. For $`\beta 1`$ charge one particles, an average of 12 photoelectrons per event were detected.

About half of the particles triggered by the instrument passed through the MWPC, where they ionized the gas. The ionization signals were amplified and detected by the pad plane along with the Cherenkov signals.

### 2.2 The time-of-flight system

The time-of-flight system consisted of two layers of plastic scintillators, one placed above the tracking system and the other below, as indicated in Figure 1. Each layer was divided into two paddles with a size of 25$`\times `$50 cm2 and a thickness of 1 cm. The material used was Bicron 401. Each paddle had two 5 cm diameter photomultiplier tubes, one at each end. The distance between the two scintillator layers was 1.2 m.

The signal from each photomultiplier was split in two parts, one was sent to an analog-to-digital converter and the other to a time-to-digital converter. In this way, the time-of-flight system provided both energy loss (dE/dX) and timing information. The scintillator signals also provided the trigger for the data acquisition system.

### 2.3 The tracking system

The tracking system consisted of a superconducting magnet and three drift chambers. The average maximum detectable rigidity (MDR) was 330 GV.

The magnet (Golden et al., 1978) consisted of a single coil of 11,161 turns of copper-clad NbTi wire. The outer diameter of the coil was 61 cm and the inner diameter 36 cm. The coil was placed in a dewar filled with liquid helium surrounded by a vacuum shell close in a second dewar filled with liquid nitrogen that reduced the rate of evaporation of liquid helium and enabled to attain a life time of about 100 hours for the superconducting magnet. The operating current was set at 120 A, producing an inhomogeneous field of approximately 4 T at the center of the coil.

The three drift chambers (Hof et al., 1994) used for the trajectory measurements were physically identical. The lateral sides of the chamber box were made from 1 cm thick epoxy-composite plates, while the open top and bottom sides were covered with 160 $`\mu `$m thick copper plated mylar windows. The inner gas volume of the box was of size $`47\times 47\times 35`$ cm3. The drift chamber had six layers, each layer containing sixteen 27.02 mm wide drift cells, for measurements in the x-directions and four layers for the y-direction. A high efficiency ($`99\%`$ for a single drift cell) and an average spatial resolution better than 100 $`\mu `$m were found. The three drift chambers provided a total of 18 position measurements in the direction of maximum bending, x direction, and 12 along the perpendicular view, y direction. Using the position information together with the map of the magnetic field, the rigidity of the particle was determined.

### 2.4 The Calorimeter

The silicon tungsten calorimeter flew in several balloon-borne experiments. The CAPRICE98 configuration also was used in the CAPRICE94 (Boezio et al., 1997) and CAPRICE97 (Kremer et al., 1999) experiments. The calorimeter (Bocciolini et al., 1996; Ricci et al., 1999) was designed to distinguish non-interacting minimum ionizing particles, hadronic cascades and electromagnetic showers.

The calorimeter consisted of eight 48 $`\times `$ 48 cm2 silicon planes interleaved with seven layers of tungsten converter, each one radiation length (X0) thick. A single plane consisted of an array of 8$`\times `$8 pair of silicon detectors. Each detector had a total area of 60$`\times `$60 mm2 and was divided into 16 strips, each of width 3.6 mm. The detectors were mounted back-to-back with perpendicular strips to give $`x`$ and $`y`$ readout. The strips of each detector were daisy-chained longitudinally to form one single strip 48 cm long. Taking into account all the material, the calorimeter had a total thickness of 7.2 X0.

## 3 Data analysis

The analysis was based on 21 hours of data for a total acquisition time of 67240 seconds under an average residual atmosphere of 5.5 g/cm2. The fractional live time during the flight was $`0.4865\pm 0.0002`$ resulting in a total live time ($`T_{live}`$) of $`32712\pm 13`$ s.

Antiprotons are a very rare component of the cosmic radiation. They must be distinguished from a large background of protons and electrons. Furthermore, products of interactions of cosmic rays in the atmosphere above the payload, such as muons and pions, are a significant background for antiproton measurements performed with balloon-borne experiments. For these reasons strict selection criteria had to be applied on the data acquired from each detector. The rigidity range of the antiproton analysis was 4 to 50 GV. The lower limit was due to the geomagnetic cut-off of the experiment, while the upper limit was based on the RICH ability to reliably identify antiprotons from other particles at maximum Cherenkov angle ($`\beta 1`$). At 50 GV the (anti)proton Cherenkov angle became less than 3 standard deviations away from the Cherenkov angle of $`\beta 1`$ particles.

### 3.1 Antiproton and proton selection

#### 3.1.1 Tracking

The primary task of the tracking system was to precisely measure the sign and absolute value of the deflection (1/rigidity) of the particle traversing the apparatus. Events with more than one track, such as products of interactions, were eliminated. For this reason a set of strict selection criteria was imposed on the quality of the fitted tracks. These criteria were based partly on experience gained during the analysis of data from a similar tracking system (Hof et al., 1996; Mitchell et al., 1996; Boezio et al., 1997):

1. At least 12 (out of 18) position measurements in the x direction (direction of maximum bending) and 8 (out of 12) in the y direction were used in the fit.
2. There should be an acceptable chi-square for the fitted track in both directions with stronger requirements on the x-direction.

In addition to these criteria, we required that the value of the deflection, as determined using only the tracking information from the top half of the spectrometer, be consistent with the value determined using only the bottom half of the spectrometer. This additional constraint was imposed on events below the RICH threshold for antiprotons. For events where Cherenkov light was detected (above 18 GV) the particles could be distinguished using the calculated mass. To remove the contamination from spillover protons in the antiproton sample, a concern at high rigidities, we required that the uncertainty in the deflection, estimated on an event-by-event basis, be less that 0.008 (GV)-1 (Golden et al., 1991). This value was chosen as a compromise between rejection power for spillover protons and efficiency of the condition in the high rigidity range. It is worth pointing out that the Cherenkov angle determined by the RICH detector provided an additional check on the deflection (see also section 3.2.5) with a comparable uncertainty at 50 GV.

#### 3.1.2 Scintillators and time-of-flight

The information of the time-of-flight system, with a time resolution of 230 ps, was used to select downward moving particles. The dE/dx information from the top scintillators was used to reject alphas and heavier particles as well as multiparticle events coming from interactions above the top scintillator. This was done by requiring the following two conditions.

1. dE/dx losses in the top scintillator less than 1.8 mip (where a mip is the most probable energy loss for a minimum ionizing particle).
2. Only one paddle hit in the top scintillator plane.

Antiprotons interacting in the calorimeter could produce backscattered particles that traverse the bottom scintillator paddles giving an additional signal. None of these cases significantly affected the performance of the tracking system and calorimeter. Therefore, no restrictions were put on the bottom scintillators.

#### 3.1.3 Calorimeter

The calorimeter was primarily used to identify electrons. The longitudinal and transverse segmentation of the calorimeter combined with the measurement of the energy lost by the particle in each silicon strip resulted in a high identification power for electromagnetic showers. In the analysis presented in this paper, the calorimeter was used to reject events with electromagnetic showers (see Bergström 2000c111Ph.D. thesis 2000, Royal Institute of Technology, is available at: http://www.particle.kth.se/group\_docs/astro/research/references.html for a description of the selection criteria), hence reducing the electron contaminations in the antiproton sample. The procedure followed was similar to the one used in the CAPRICE94 antiproton analysis (Boezio et al., 1997). The selection was designed to reject electrons while keeping as large an antiproton fraction as possible.

Figure 2 illustrates the calorimeter performance and shows a schematic view of a 5 GV electron in the CAPRICE98 apparatus. In the figure there are a left- and a right-view, symbolizing respectively the x and y views of the CAPRICE98 apparatus. At the top is the RICH detector and a rotated view of the signals in the pad plane of the multiwire proportional chamber is showed in the square frame in the center of the figure. The ionization cluster of pads can be seen well separated from the Cherenkov ring typical of a $`\beta 1`$ particle. The three central boxes are the drift chambers of the tracking system. At the bottom there is a rectangular frame that symbolizes the calorimeter. The line that is drawn through all detectors represents the fitted track of the particle. Along the line in the drift chambers there are small circles drawn around each wire that gave a signal. The size of the circle is proportional to the calculated drift time for the electrons at that wire. The calorimeter shows the electromagnetic shower induced by an electron.

The electromagnetic shower in the calorimeter of Figure 2 is clearly distinguishable from the hadronic shower produced by an interacting antiproton candidate shown in Figure 3 and from the non-interacting pattern of another antiproton candidate shown in Figure 4.

Out of the 31 antiproton events surviving all antiproton selections 8 were found to have interacted in the calorimeter. This is in agreement with the simulated expectation of $`9.7\pm 1.7`$ antiproton interactions in the calorimeter.

#### 3.1.4 RICH

The RICH was used to measure the Cherenkov angle of the particles and hence their velocities. Below the threshold rigidity for antiprotons to produce Cherenkov light in the gas (about 18 GV) the detector was used as a threshold device to veto lighter particles, while above it the Cherenkov angle was reconstructed. Below 25 GV the fluctuations in the number of detected photoelectrons were quite large due to variations in threshold rigidity caused by pressure variation in the radiating gas and because of large Poisson fluctuations as the average number of photoelectrons detected at this rigidity was 6. Therefore antiprotons in the rigidity range 4 to 25 GV were selected if the event did not produce a Cherenkov signal or when the reconstructed Cherenkov angle was consistent with that of an antiproton with the measured rigidity.

We show in Figure 5 the separation between antiprotons and lighter particles, after applying tracking, ToF and calorimeter selection criteria, by means of the RICH information in the rigidity range from 4 to 18 GV. The events in Figure 5 are plotted as a function of the logarithm in base 10 of the number of pads ($`n`$) used for the reconstruction of the Cherenkov angle (Bergström 1999222licentiat thesis 1999, Royal Institute of Technology, is available at: http://www.particle.kth.se/group\_docs/astro/research/references.html; Bergström et al. 2000b) plus 1. It can be noted from the top panel that protons have low values of $`n`$ (note the logarithmic scale on the Y axis), 93% being at zero, while on the negative side a corresponding peak can be seen indicating the antiprotons. It is evident from the figure that faster particles (mostly muons and pions) at higher $`n`$ can clearly be separated from the antiprotons. In the analysis, (anti)protons below 25 GV were selected by requiring a value of $`n`$ equal to 0. In this way 26 antiprotons were identified. These events had rigidities between 4 and 17 GV. Note that events were also selected if the Cherenkov angle could be constructed and 3 events were identified below 25 GV.

It is worth noting that when using the RICH detector as a threshold device the tracking system gave information about where in the pad plane the Cherenkov light should have been detected, thereby greatly enhancing the detector immunity to noise. The number of noisy pads per event was, on the average, less than one out of 4096 channels/pads. This allowed a stringent selection to be applied for antiprotons where no signals should appear in the predicted area, and yet maintaining a high identification efficiency. All 4096 channels were working during the entire flight.

Figure 3 shows the schematic view of a 6.4 GV selected antiproton. At this rigidity the antiproton does not produce light and, in fact, only the signal from the ionization due to the crossing particle is detected in the pad plane.

Since the RICH detector played a crucial role in the rejection of background muon and pion events, stringent selection criteria were applied to the RICH data for the events above the antiproton threshold. The Cherenkov angle resolution for protons was determined using a large sample of protons selected using the calorimeter and scintillators. The resolution varied from 11 mrad at threshold to about 1.3 mrad for fully relativistic protons Bergström et al. (2001). Figure 6 shows the measured Cherenkov angle as a function of rigidity. The events were selected from the flight data after applying tracking and ToF selection criteria. The solid, dotted and dashed lines indicate the calculated Cherenkov angle for muons, kaons and (anti)protons respectively. Five antiprotons are clearly identified and they are shown with black boxes and with one standard deviation error bars for both the rigidity and Cherenkov angle measurements.

In conclusion the conditions on the RICH detector information used for the antiproton selection were:

1. The center, extrapolated from the tracking information, of the Cherenkov ring was required to be contained in the pad plane.
2. Multiple charged tracks traversing the MWPC were rejected by requiring that there was only one cluster of pads with a high signal, typical of ionization from a charged particle in the location indicated by the tracking system.
3. If the particle crossed the MWPC (46% of the events), a good agreement between the particle’s impact position as determined by the RICH and the tracking system was required. The difference in x and y should be less than 3 standard deviations (rigidity dependent), typically $`<`$ 5 mm.
4. Conditions on the Cherenkov angle for events:
(i) Between 4 and 25 GV, there was no signal due to Cherenkov light.
(ii) However, the following criteria were met for all rigidities above the antiproton threshold (calculated on an event by event base according to the measured gas pressure (Bergström et al., 2001)):

1. A rigidity dependent condition on the number of pads used for the reconstruction of the Cherenkov angle was applied. The condition required more than 5 pads at the antiproton threshold increasing to 20 above 35 GV.
    2. The reconstructed Cherenkov angle should not deviate by more than 3 standard deviations below and 2 standard deviations above from the expected Cherenkov angle for (anti)protons.
    3. To suppress the background from lighter particles, the reconstructed Cherenkov angle was required to be more than 4 mrad (3 standard deviations for $`\beta 1`$ particles) away from the expected Cherenkov angle for pions (about 53 mrad above 18 GV).

Hence between the antiproton threshold (about 18 GV) and 25 GV, antiprotons were selected either if the event did not produce a Cherenkov signal or if the reconstructed Cherenkov angle was consistent with that of an antiproton with the measured rigidity.

Figure 4 shows the schematic view of one of the selected antiprotons. The clean Cherenkov ring in the RICH pad plane is well separated from the ionization cluster of pads, hence permitting the mass associated with the event to be reconstructed.

### 3.2 Contamination

The contamination due to e-, $`\mu ^{}`$, $`\pi ^{}`$ and spillover protons in the antiproton sample was studied carefully using simulations and experimental data taken during the flight and on the ground before the flight.

#### 3.2.1 Albedo particles

Albedo particles were rejected using the time-of-flight information. With a time-of-flight between the top and bottom scintillators of more than 4 ns, the 0.23 ns resolution ensured a negligible background of upgoing particles.

#### 3.2.2 Electron contamination

The calorimeter performance was studied primarily with simulations. Simulation studies showed that electron contamination in the calorimeter selection was of ($`0.6\pm 0.2`$)% independent of rigidity in the interval from 4 to 50 GV. This value was cross-checked by studying electrons selected using a condition on the total number of strips hit in the calorimeter. It is worth noting that only the calorimeter was able to separate electrons from muons above about 5 GV.

The electron contamination in the RICH selection was studied using a sample of 495 e- in the interval 4 to 50 GV (478 between 4 and 18 GV), selected using the calorimeter. Of the 495 events one was selected as antiproton resulting in an electron contamination of ($`0.20_{0.17}^{+0.47}`$)% in the RICH selection for antiprotons.

To estimate the electron background, we considered the events with negative curvature between 4 and 50 GV. After imposing the tracking and ToF selection criteria to this sample, 1031 negative events were left. As a worst case, we assumed that all of these events were electrons. Applying the rejecting power of the RICH and calorimeter to this sample resulted in an electron contamination of less than 0.1 event over the entire range of the antiproton measurement. Hence, the electron contamination in the antiproton sample can be assumed negligible.

#### 3.2.3 Muon contamination

The calorimeter cannot separate muons from non-interacting antiprotons. The antiproton identification in a muon background was performed by the RICH. The muon contamination in the RICH selection was studied using a sample of muons collected during a ground data run prior to the launch. The fraction of muons surviving the antiproton RICH selection was ($`0.36_{0.10}^{+0.13}`$)% between 4 and 8 GV, ($`0.30_{0.12}^{+0.18}`$)% between 8 and 18 GV, ($`0.76_{0.36}^{+0.60}`$)% between 18 and 30 GV and ($`1.7_{0.8}^{+1.4}`$)% between 30 and 50 GV. Defining negative muons all the events surviving the tracking, ToF and calorimeter antiproton selection criteria, 319, 108, 24 and 11 muons were selected from the flight data between 4 and 8 GV, 8 and 18 GV, 18 and 30 GV and 30 and 50 GV, respectively. Multiplying these numbers by the surviving fractions found above and taking into account the presence of antiprotons in the sample, the muon contamination in the antiproton sample was estimated to be $`1.1_{0.3}^{+0.4}`$ in the first rigidity bin, $`0.28_{0.11}^{+0.17}`$ in the second, $`0.15_{0.07}^{+0.12}`$ in the third and $`0.10_{0.08}^{+0.15}`$ in the forth. This contamination was later subtracted from the antiproton signal and is shown in parenthesis in Table 1.

#### 3.2.4 Meson contamination

Pions started to produce light in the gas-RICH above about 3 GV. Since the RICH is a $`\beta `$-detector, the pion contamination was studied by scaling the muon sample. The contamination was similar to that of muons except in the first bin where it was ($`0.6\pm 0.1`$)%. However, pions were a small component compared to muons in the flight data. Between 3 and 4 GV pions could be clearly identified by the gas-RICH and it was found that pions were ($`12\pm 2`$)% of the muons. Assuming that a similar value is valid for the bin 4 to 8 GV, the resulting pion plus muon contamination in the antiproton sample in this bin was $`1.17_{0.27}^{+0.35}`$.

Theoretical calculations (Stephens, 1981) of the kaon spectrum in the atmosphere indicate that the kaon to antiproton ratio is about 2%. Hence, they are a negligible contamination of the antiproton sample. Furthermore, above about 9 GV kaons are suppressed by the RICH selection. It is worth pointing out that, as can be seen in Figure 6 (dotted line) , no negative kaons were identified between 9 and 50 GV.

#### 3.2.5 Spillover proton contamination

Spillover protons can represent a non negligible contamination in the antiproton sample above 20 GV.

To obtain the shape of the spillover proton distribution a Monte Carlo approach was used. We started with an input power law spectrum in rigidity that was transformed in deflection and smeared with values randomly picked from the resolution function. The spectral index for the power law was obtained from the proton spectrum measured by CAPRICE98 above 20 GV and was found to be: $`2.74\pm 0.02`$.

The resolution function was obtained from a measurement of straight tracks taken at ground prior to the flight without the magnetic field. These were analyzed as if they were high rigidity events with magnet on. In this case the resolution function is simply the deflection distribution. To construct this distribution events that suffered multiple Coulomb scattering had to be eliminated since they could enlarge the distribution (see Menn et al., 2000). In fact, the multiple scattering affected events essentially at low energy while spillover events were due to high energy protons. This was done by requiring a Cherenkov light signal along with high value of the measured Cherenkov angle in the RICH so as to select close to fully relativistic events (mostly muons and electrons). Then electrons were eliminated using the calorimeter. The large majority of the selected events were particles (muons) with rigidity greater than 5 GV.

The simulated proton spillover distribution was normalized with the experimental deflection distribution between -0.02 and 0 (GV)-1 and used to obtain an estimation of the experimental proton spillover distribution. Figure 7 shows the experimental distribution with the estimated proton spillover contribution (solid line). The distribution includes spillover protons, muons and antiprotons. The events for the deflection distribution were selected from the whole flight data set with the complete antiproton selection (tracking, ToF, RICH and calorimeter conditions) except for the condition on the reconstructed Cherenkov angle. At the smallest deflections (high rigidities) the dominant component was the spillover protons.

Finally, the contamination of spillover protons was determined from this calculation by integrating the experimental spillover distribution over the proper deflection range. Between 0 and 30 GV the estimated spillover contamination was less than 0.002 events and between 0 and 50 GV it was about 0.39. However, this was not the real spillover contamination in the antiproton sample since antiprotons above the threshold of the RICH detector were selected with the additional condition on the measured Cherenkov angle. For the rigidity range from the threshold of the RICH detector up to 50 GV the spillover proton contamination was obtained by calculating the probability that the 0.39 spillover protons previously estimated were selected with the Cherenkov angle condition. This probability was obtained using negative events with deflection smaller than 0.01 GV-1 (i.e. essentially spillover protons). Using this process, we derived that the contamination of spillover protons was 0.04 events below 50 GV.

### 3.3 Efficiency

The antiproton selection efficiencies were studied using a large sample of experimental protons from flight data set. It was assumed that protons and antiprotons had the same efficiencies in the RICH, scintillators, and the tracking system. However, in the calorimeter the efficiencies were not assumed to be the same because of the difference in the interaction cross sections. Hence, the calorimeter efficiency for protons and antiprotons was studied using simulations.

The tracking efficiency was studied with two independent methods (see Bergström et al. (2000)). Both were based on track reconstruction codes that were independent of the drift chamber tracking system and used information from the other available detectors. The first method used the position of the ionization cluster in the MWPC of the RICH detector, information from the two ToF scintillation detectors (only in the x direction) and from the electromagnetic calorimeter to reconstruct the track of the particle traversing the detector system. This combination had an estimated MDR of 4.5 GV. The second method used the RICH to determine the rigidity from the velocity derived from the Cherenkov angle measurement with the help of an extrapolated straight track from the calorimeter.

The two methods were tested with muons from ground data and similar results were obtained between 0.2 and 10 GV. With flight data, the first method sampled the tracking efficiency of protons below 10 GV while the second above 18 GV (because of the gas-RICH threshold). In this case the efficiencies differed by $`6\%`$. From previous experience with a similar tracking system (Boezio et al., 1999) the proton tracking efficiency was expected to reach a plateau above 2 GV. This seems in disagreement with what found here. However, the first method could be biased by contamination of secondary low energy protons that would not affect the second method (for more details see Bergström et al. (2000)). Thus, the efficiency of the tracking selection was obtained by the second method for the rigidity range from 4 to 50 GV (dotted line in Figure 8). The difference found between the proton tracking efficiency with the two methods was considered a systematic uncertainty and a 6% systematic uncertainty was included in the flux calculation for rigidities between 4 and 20 GV.

Secondary particles backscattering from the calorimeter could result in an inefficiency for the tracking selection. We studied this effect selecting from the flight data events that interacted in the calorimeter producing large signals (above 1.8 mip) in the bottom scintillator. For these events we applied the tracking selection. Accounting for the relative abundance of these events, we found an overall decrease in the tracking efficiency of about 0.4%.

A Monte Carlo simulation based on the CERN GEANT/FLUKA-3.21 code (Brun et al., 1994) was used to study the calorimeter selection efficiencies. Results from similar simulations were found to be in good agreement with test beam data and experimental results from previous balloon flights (Bocciolini et al. 1993; Boezio 1998333Ph.D. thesis 1998, Royal Institute of Technology, is available at: http://www.particle.kth.se/group\_docs/astro/research/references.html). It was found that the efficiency for the selection of antiprotons and protons using simulated data were in agreement, within the statistical uncertainties. The calorimeter efficiency for protons selected from experimental data was in reasonable agreement with previous simulations. Differences (about 2%) found were probably due to the less accurate simulations of hadronic showers compared to electromagnetic showers. Since simulations indicated that protons’ and antiprotons’ calorimeter selections had the same efficiency, the efficiency of the experimental proton calorimeter selection was also used for the antiprotons (dashed line in Figure 8).

The RICH efficiency was rigidity dependent and is shown as solid line in Figure 8. As expected the efficiency was constant below the gas-RICH threshold but it started decreasing at 14 GV due to above threshold protons that spilled down to lower rigidities because of the finite resolution of the tracking system. Above 25 GV the antiprotons were selected only with conditions on the Cherenkov angle . The decrease of the RICH efficiency above 30 GV was caused by the requirement that the Cherenkov angle should be more than 4 mrad away from the expected Cherenkov angle for pions.

Since the detector efficiencies varied with rigidity, mean efficiencies had to be calculated. This was done by weighting the efficiencies in each bin with the proton spectrum measured in this experiment (with a larger number of bins and no RICH selection) and with an antiproton spectrum given by interstellar secondary calculation by L. Bergström & P. Ullio (1999, private communication) which included the effects if the geomagnetic cut-off.

### 3.4 Geometrical factor

The acceptance of the instrument allowed for particles with a range of zenith angles to be measured. The maximum angle was 14 degrees and the mean of the distribution was at 8 degrees.

The geometrical factor was obtained with Monte Carlo techniques (Sullivan, 1971) and the same track-fitting algorithm used in this analysis to trace the particles through the spectrometer. The geometrical factor ($`G`$) was found to be constant at $`155.0\pm 1.1`$ cm2 sr in the rigidity range from 4 to 50 GV.

### 3.5 Payload and atmospheric corrections

#### 3.5.1 Payload corrections

To determine the number of particles at the top of the payload, the losses and the production of particles due to interaction in the material of the payload had to be considered. To reach the tracking system of the spectrometer, the particles had to go through first the aluminum shell of the payload, the RICH detector and then the top scintillator of the time-of-flight system. It was assumed that all particles that interacted above the tracking system were rejected by the selection criteria. The probability of an interaction in the material of the drift chambers, that would not be rejected by the tracking system conditions, was considered as being negligible. The data were corrected for these losses with multiplicative factors, using the expression for the interaction mean free path for the different materials in the detectors given by Stephens (1997). This gave correction factors of 1.132, 1.118, 1.109 and 1.104 for antiprotons and of 1.081, 1.082, 1.084 and 1.086 for protons in the four rigidity intervals. The corrected number of antiprotons and protons at the top of the payload are given in Tables 1 and 3.

#### 3.5.2 Atmospheric corrections

For the production of secondary antiprotons and protons in the atmosphere, we used the calculation by Papini, Grimani, & Stephens (1996) for the protons and the calculation by Stephens (1997) for the antiprotons. To determine the secondary spectra we compared our measured proton spectrum propagated to the top of the atmosphere with the spectrum used for the secondary calculations at solar minimum (Papini, Grimani, & Stephens, 1996; Stephens, 1997) from which a normalization factor was derived. The resulting secondary fluxes were normalized with the geometrical factor and live time of the experiment, and subtracted from the corrected numbers using a mean residual atmosphere of 5.5 g/cm2. The number of atmospheric antiprotons and protons are given in Tables 1 and 3.

The correction for losses in the atmosphere was carried out using a method analogous to the instrument correction. This gave correction factors of 1.1, 1.089, 1.082 and 1.079 for antiprotons and of 1.061, 1.063, 1.064 and 1.065 for protons in the four rigidity intervals.

### 3.6 Geomagnetic transmission correction

To be able to get the fluxes at the top of the atmosphere, the transmission of the particles through the earth’s magnetic field had to be taken into account. During the flight the position of the payload changed between 34.3 and 35.5 North Latitude and between 104.1 and 110.1 West Longitude. This correspond to an average value of the vertical cut-off rigidity of about 4.3 GV (Shea & Smart, 1983). However, this cut-off is not a sharp value below which all particles are deflected and cannot reach the apparatus and above which all particles arrive. In fact, around the geomagnetic cut-off the particles are partially transmitted through the earth magnetic field. Furthermore, the penumbral bands define forbidden bands of rigidity, which vary with arrival direction, time and geographic location. In this analysis all these effects are represented by a transmission function which was derived by the experimental data.

We found that the CAPRICE94 (Boezio et al., 1999) and CAPRICE98 proton spectra above about 10 GeV are nearly identical in shape and the absolute fluxes differ by $`7\%`$ in good agreement considering both the statistics and systematic errors, which in Boezio et al. (1999) were estimated to be of the order of 10%. Moreover, the solar modulation during the two balloon flights was also very similar. The values from the neutron monitor counter CLIMAX444National Science Foundation grant ATM-9912341, http://ulysses.uchicago.edu/NeutronMonitor/Misc/neutron2.html (Simpson & Pyle, 1996) were 415600 counts/hour and 417000 counts/hour at the time of the CAPRICE94 and CAPRICE98 flights, respectively. However, the CAPRICE94 experiment took place in North Canada at an average geomagnetic cut-off of about 0.5 GV. Hence, above 1 GV the effects of the geomagnetic field on the CAPRICE94 proton spectrum were negligible. Consequently, the transmission function was defined as the ratio between the experimental CAPRICE98 and CAPRICE94 proton fluxes.

The correction factors ($`TF`$) for the geomagnetic effect were derived weighting the transmission function with the proton and antiproton spectra as done for the efficiencies. The resulting transmission values differ from 1 only in first bin (4 to 8 GV) where they are $`0.84\pm 0.06`$ and $`0.81\pm 0.07`$ for antiproton and protons, respectively.

## 4 Results

### 4.1 Antiproton flux at the Top of the Atmosphere

Given the number of events ($`N_{\overline{\mathrm{p}}\text{ }}^{TOA}`$) selected with the antiproton criteria and corrected for selection efficiencies, losses in the payload and in the atmosphere and atmospheric secondaries, we obtained the antiproton fluxes at the top of the atmosphere from the relation,

$$\text{Flux(E)}=\frac{1}{T_{live}\times G\times \mathrm{\Delta }E\times TF}\times N_{\overline{\mathrm{p}}\text{ }}^{TOA}(E),$$

where $`\mathrm{\Delta }E`$ is the energy bin corrected for ionization losses to the top of the atmosphere and $`E`$ the kinetic energy. The resulting antiproton flux is given in Table 2. The total errors include both statistical and systematic errors. The mean energies of the bins are given according to Lafferty & Wyatt (1995).

Figure 9 shows the antiproton flux measured by this experiment together with other experimental data (Buffington, Schindler, & Pennypacker, 1981; Mitchell et al., 1996; Boezio et al., 1997; Basini et al., 1999; Orito et al., 2000). The two solid lines show the upper and lower limit of a calculated flux of interstellar secondary antiprotons (Simon, Molnar, & Roesler, 1998) using a recently measured proton and helium spectra (Menn et al., 1997) and a reacceleration model which allows energy-changing processes caused by the non-annihilation process and by elastic scattering. The dashed line shows the interstellar secondary antiproton flux calculated by L. Bergström & P. Ullio (1999, private communication). This calculation assumed a diffusion model of propagation with an isotropic diffusion coefficient and no reacceleration. It used the interstellar proton spectrum measured by the CAPRICE94 experiment (Boezio et al., 1999). The dotted line shows the primary antiproton flux by Ullio (1999), which included a Minimal Supersymmetric Standard Model with a contribution from an assumed Higgsino-like neutralino, with a mass of 964 GeV. The theoretical fluxes, but not the experimental values of the other experiments, were corrected for the solar modulation conditions corresponding to the CAPRICE98 flight using a spherically symmetric model (Gleeson & Axford, 1968) with a solar modulation parameter of $`\mathrm{\Phi }=600`$ MV.

### 4.2 Antiproton to Proton ratio

To obtain the antiproton to proton ratio at the top of the atmosphere, we corrected the number of selected antiprotons and protons for the production and loss of particles in the residual atmosphere above the apparatus as well as in the instrument itself. However, for obtaining antiproton to proton ratio, the selection efficiencies, which were considered to be the same for antiprotons and protons, were excluded from the calculation in order to reduce the errors. The resultant ratios are presented in Table 3.

Figure 10 shows the antiproton to proton ratio measured by CAPRICE98 along with other experimental data (Buffington, Schindler, & Pennypacker, 1981; Golden et al., 1984; Bogomolov et al., 1987, 1990; Salomon et al., 1990; Stochaj et al., 1990; Mitchell et al., 1996; Boezio et al., 1997; Basini et al., 1999; Orito et al., 2000) and with theoretical calculations. The two solid lines show the upper and lower limit of a calculated flux of interstellar antiprotons by Simon, Molnar, & Roesler (1998) assuming a pure secondary production during the propagation of cosmic rays in the Galaxy . The dashed line shows a similar calculation by L. Bergström and P. Ullio (1999, private communication). It is worth noting that Simon, Molnar, & Roesler (1998) used the primary spectra measured by Menn et al. (1997) and that L. Bergström and P. Ullio (1999, private communication) used the interstellar proton spectrum measured by Boezio et al. (1999).

The antiproton to proton ratio values presented here are in perfect agreement with the preliminary CAPRICE98 results published by Bergström et al. (2000). A small variation, of about a third of standard deviation, is found only for the highest energy bin essentially due to a better understanding of the proton spillover contamination.

## 5 Systematic uncertainties

Systematic errors originating from the determination of the detector efficiencies have already been discussed and they have been included in the data in the figures and in the tables.

Another possible systematic error is related to the efficiency of the trigger system. The trigger efficiency was studied during the pre-flight preparations with a system measuring the coincidence of two scintillators placed above and below the top and bottom time-of-flight scintillators. The efficiency was found to be close to 100% with an uncertainty of about 2%. The performance of the trigger system during the flight was also studied comparing the experimental spatial distribution of triggers with the distribution given by the same simulations as were used for the geometrical factor calculation, and an excellent agreement was found. Hence, the trigger efficiency could be assumed to be 100%, with a possible systematic uncertainty of less than 2%.

The method for calculating the geometrical factor used in this work was compared with two other techniques in the CAPRICE94 analysis (Boezio et al., 1999), and it was found to be in agreement within 2%, above 0.5 GV. Considering the similar geometrical configuration of CAPRICE98, it can be concluded that the uncertainty on values of the geometrical factor was about 2%.

Systematic errors due to the uncertainty on the secondary production of antiprotons in the atmosphere were estimated comparing the results from the calculation by Stephens (1997), used in this work, with the independent calculation by Pfeifer, Roesler, & Simon (1996). It was found that the two calculations differ by $`22\%`$ at 5 GeV decreasing to $`2\%`$ above 8 GeV, with the calculation by Pfeifer, Roesler, & Simon (1996) being the lower. This introduces an estimated uncertainty in the antiproton fluxes extrapolated to the top of the atmosphere, that is $`6\%`$ at 5 GeV decreasing to less than 1% above 8 GeV.

The atmospheric secondaries were also affected by the uncertainty in the residual atmosphere above the gondola. This was measured to be 5.5 g/cm2 by a pressure sensor owned and calibrated by the CAPRICE collaboration. The pressure was also measured by a detector owned and calibrated by the National Scientific Balloon Facility (NSBF). The NSBF pressure data were about 15% lower at float than the ones measured by our sensor. This results in an uncertainty on the antiproton fluxes of about 7% between 3 and 20 GeV decreasing to about 2% above 30 GeV.

The numbers of particles measured at the spectrometer were corrected for losses in the spectrometer and the atmosphere. Assuming a 10% uncertainty on the cross sections used in these calculations results in a systematic error on the antiproton fluxes of $`2\%`$. An additional uncertainty of 1% should be considered due to the uncertainty on the atmospheric depths and the consequent effect on the losses in the atmosphere.

Assuming that the systematic errors were uncorrelated and, hence, could be quadratically summed, the resulting measurement of the antiproton flux includes systematic uncertainties that were energy dependent decreasing from $`10\%`$ at 5 GeV to $`8\%`$ above 8 GeV and to $`4\%`$ above 30 GeV.

Systematic uncertainties due to the tracking system, caused by, for example, an offset in the deflection measurements, were analyzed using the RICH detector. The RICH high Lorentz threshold for protons permitted to study several features of the tracking system up to a rigidity of about 100 GV (unattainable by previous cosmic rays experiments).

Figure 11(a) shows the Cherenkov angle resolution for protons obtained from flight data ($``$) as a function of $`\beta `$, derived from the rigidity measured by the tracking system assuming the proton mass. The solid line indicates the measured resolution for muons (Bergström et al., 2001). The difference between the two resolutions was due to the finite resolution of the tracking system. In fact, this introduced an additional spread in the Cherenkov angle distribution, when it was binned according to the $`\beta `$ measured by the tracking system. The effect of this was more important for protons than for muons, since the rigidity (deflection) of the protons was nearly 10 times greater (smaller) than that of the muons, at the same $`\beta `$.

The effect of the tracking resolution of the binned Cherenkov angle distribution was obtained by simulating a large number of protons according to a power law spectrum in rigidity with spectral index of -2.74. The corresponding Cherenkov angle was then calculated and smeared with a gaussian distribution with a standard deviation given by the measured resolution for muons (Fig. 11(a), solid line). Then the rigidity of each simulated proton was transformed to a corresponding deflection and smeared with values randomly picked from the tracking resolution function. The resulting deflection was then used to derive the velocity, which was used for the binning of the Cherenkov angles of these simulated events, similarly to the experimental case. The resulting Cherenkov angle resolution is shown as $`\mathrm{}`$ in Figure 11 (a). A good agreement was found between the measured and simulated resolutions. This was a strong indication that the tracking resolution function used in this work describes the tracking uncertainties with high precision and, consequently, could be used for determining the spillover proton contamination in the antiproton sample.

Figure 11(b) shows the measured mean Cherenkov angle ($``$) as a function of $`\beta `$ from the tracking system along with the simulated ($`\mathrm{}`$) one. The same simulation was used as for the Cherenkov angle resolution. Also in this case, an excellent agreement was found. Furthermore, the comparison between the two sets of mean Cherenkov angle permitted to limit the possibility of an offset in the measured deflection, due to effects such as: wrong alignment of the drift chambers, positioning of the center of the magnet, mapping of the magnetic field, etc. An offset was introduced in the simulated deflection, and it was varied over a wide range of possible values. It was found that, if an offset existed, at a 95% confidence level it was not larger than 0.001 (GV)-1, which was significantly smaller than 0.003 (GV)-1, corresponding to the MDR of the experiment (330 GV).

It is worth pointing out that the simulation was tested also with different spectral index such as: -2.6 and -2.8, and no significant variation from the case presented here was found.

## 6 Conclusion

The antiproton flux and the antiproton to proton ratio have been determined in the energy region from 3 to 49 GeV by the CAPRICE98 experiment. This is the first time that the antiproton flux has been measured up to such high energies and over such a wide range in energy. Between 3 and 20 GeV our antiproton fluxes are consistent with the measurement by the MASS91 experiment (Basini et al., 1999). Both of these results, within the experimental errors, are also in agreement with the theoretical predictions by Simon, Molnar, & Roesler (1998) and L. Bergström and P. Ullio 1999 (private communication) which assume a purely secondary origin of the cosmic-ray antiprotons. However, a primary component cannot be excluded and the shape of the measured antiproton flux could indicate a presence of primary antiprotons. In fact, in the CAPRICE98 analysis we observed two antiproton events, with the highest energy antiproton measured at a kinetic energy of 43 GeV, between 29 and 49 GeV, while the expected number from a pure secondary origin is only 0.2 to 0.4 events, including muon contamination from the atmosphere; the lower and upper values correspond to the two extreme secondary curves of Figure 9. It is essential to improve the statistics on antiproton measurements in this high-energy region In fact, the energy region studied here permits to search for specific signatures of neutralino-induced antiproton fluxes which are not attainable in lower energy regions. Furthermore, in this energy range nearly all calculations of interstellar secondary antiprotons in the literature are consistent with each other (Ullio, 1999). These data are also substantially free of uncertainties due to solar modulation effects.

As a continuation of its ballooning activity, the WIZARD collaboration has developed a cosmic-ray experiment, named PAMELA (Adriani et al., 1999), that will be launched in quasi-polar orbit on board of a Russian satellite at the beginning of 2003. PAMELA is based on a magnetic spectrometer with an MDR exceeding 800 GV and will determine the antiproton spectrum from 80 MeV, with a few orders of magnitude better statistics at energies above 5 GeV than the existing ones, free of the atmospheric background. Another cosmic-ray experiment, AMS (Ahlen et al., 1994), dedicated to the search for antinuclei, will be installed sometime later on the International Space Station. As it will have a similar MDR and will be based on the same set of detectors as PAMELA, it also will give further results on the study of high energy antiprotons. Its much larger acceptance also will allow improvement to the statistics in the same high energy range that will be explored by PAMELA.

This work was supported by NASA Grant NAGW-110, The Istituto Nazionale di Fisica Nucleare, Italy, the Agenzia Spaziale Italiana, DARA and DFG in Germany, the Swedish National Space Board and the Knut and Alice Wallenberg foundation. The Swedish-French group thanks the EC SCIENCE programme for support. We wish to thank the National Scientific Balloon Facility and the NSBF launch crew that served in Fort Sumner. We would also like to acknowledge the essential support given by the Gas Work Group (EST/SM/SF) and the Thin Films & Section (EP/TA1/TF) at CERN, the LEPSI and CRN-Strasbourg and the technical staff of NMSU and of INFN.


format: text

text:
# Robust Patterns in Food Web Structure

## Abstract

We analyze the properties of seven community food webs from a variety of environments—including freshwater, marine-freshwater interfaces and terrestrial environments. We uncover quantitative unifying patterns that describe the properties of the diverse trophic webs considered and suggest that statistical physics concepts such as scaling and universality may be useful in the description of ecosystems. Specifically, we find that several quantities characterizing these diverse food webs obey functional forms that are universal across the different environments considered. The empirical results are in remarkable agreement with the analytical solution of a recently proposed model for food webs.

In natural ecosystems species are connected through trophic relationships Martinez00 ; McKann98 ; Paine92  defining intricate networks Strogatz01 ; Watts98 ; Amaral00 ; Newman01 , the so-called food webs. Understanding the structure and mechanisms underlying the formation of these webs is of great importance in ecology ecology . For this reason, much research has been done in constructing empirical webs and uncovering unifying patterns describing their structure Lawton88 ; ecology . However, in the last decade the construction of larger and more complete food webs clearly indicated that the previously reported unifying patterns do not hold for the new webs Hall ; Hall2 . Indeed, the complexity of the new webs has rendered quite difficult the challenge to obtain quantitative patterns that substitute the old ones.

Here, we analyze the properties of seven detailed community food webs from a variety of environments—including freshwater habitats, marine-freshwater interfaces and terrestrial environments. Remarkably, we uncover quantitative unifying patterns that describe the properties of the diverse trophic webs considered and capture the random and non-random aspects of their structure. Specifically, we find that several quantities—such as the distributions of number of prey, number of predators, and number of trophic links—characterizing these diverse food webs obey robust functional forms that depend on a single parameter, the linkage density $`z`$.

In our analysis, we use results obtained for complex networks Newman01  and for a recent model of food web formation, the “niche model” of Ref. Martinez00 . We first describe the theoretical model and its predictions: An ecosystem with $`S`$ species and $`L`$ trophic interactions between these species, defines a network with $`S`$ nodes and $`L`$ directed links. In the niche model, one first randomly assigns species $`i=1,\mathrm{},S`$ to “trophic niches” $`n_i`$ which are mapped into the interval . A species $`i`$ is characterized by its niche parameter $`n_i`$ and by its list of prey. Prey are chosen according to the following procedure: species $`i`$ preys on the species $`j`$ with niche parameters $`n_j`$ inside a segment of length $`xn_i`$ centered in a position chosen randomly inside the interval $`[xn_i/2,n_i]`$. Here, $`0x1`$ is a random variable with probability density function $`p_x(x)=b\left(1x\right)^{\left(b1\right)}`$ Martinez00 . The values of the parameters $`b(S^2/2L1)`$ and $`S`$ determine the linkage density $`zL/S`$ of the food web, and the directed connectance $`L/S^2`$, which is a measure of the fraction of the actual number of trophic links as compared to the maximum possible number Martinez00 .

In the limit of large web sizes ($`S1`$) and small connectances ($`L/S^21`$), one can derive analytical expressions for the distribution of number of prey $`k`$ Camacho01 . We consider the cumulative distribution $`P_{\mathrm{prey}}(k)=_{k^{}k}p_{\mathrm{prey}}(k^{})`$ because it is less noisy than the probability function $`p_{\mathrm{prey}}(k)`$. We obtain

$$P_{\mathrm{prey}}(k)=\mathrm{exp}\left(\frac{k}{2z}\right)\frac{k}{2z}E_1\left(\frac{k}{2z}\right),$$
(1)
where $`E_1(x)`$ is the so-called exponential-integral function numrec . Equation (1) predicts that the distribution of number of prey decays exponentially for large $`k`$.

Also in the limit of large web sizes and small connectances, one can derive analytical expressions for the distribution of number of predators $`m`$ Camacho01 . We obtain,

$$P_{\mathrm{pred}}(m)=\frac{1}{2z}\underset{m^{}=m}{\overset{\mathrm{}}{}}\gamma (m^{}+1,2z),$$
(2)
where $`\gamma (m+1,z)`$ is the so-called “incomplete gamma function” numrec . To gain intuition about the functional form (2), note that $`p_{\mathrm{pred}}(m)`$ is approximately a step function: It is constant for $`m<z`$, and then it decays with a Gaussian tail for $`m2z`$ Camacho01 . It follows then that the cumulative distribution $`P_{\mathrm{pred}}(m)`$ decreases linearly as $`1m/z`$ for $`m<z`$ and decays as the error function numrec  for $`m2z`$.

Next, we analyze the empirical data for seven food webs with 25 to 92 trophic species. These webs have linkage densities $`2.2<z<10.8`$, and connectances in the interval 0.06–0.31 Martinez00 . We first investigate the distributions of number of prey and number of predators. Figures 1a,b compare the cumulative distributions of the number of prey and number of predators for species in the St. Martin Island web smartin  with our analytical predictions, and suggest that these distributions are well approximated by Eqs. (1)–(2) without any free parameters for fitting. Equations (1)–(2) and the results of Fig. 1 suggest the possibility that $`P_{\mathrm{prey}}`$ and $`P_{\mathrm{pred}}`$ obey universal functional forms that depend only on $`z`$.

Indeed, Eq. (1) predicts that $`P_{\mathrm{prey}}(k)`$ depends only on $`k/2z`$. So, we plot in Figs. 2a,c the cumulative distributions $`P_{\mathrm{prey}}(k)`$ versus the scaled variable $`k/2z`$ for the food webs and find that the data collapse onto a single curve, supporting the possibility that $`P_{\mathrm{prey}}`$ obeys a universal functional form ks-test1 .

The scaling of $`P_{\mathrm{pred}}(m)`$ is not as straightforward. Equation (2) indicates that “true” scaling holds only for $`m/2z<1/2`$, while for larger values of $`m/2z`$ there is a Gaussian decay of the probability function with an explicit dependence on $`z`$. However, the decay for $`m>2z`$ is quite fast and, to first approximation, not very relevant. Thus, we plot $`P_{\mathrm{pred}}(m)`$ versus the scaled variable $`m/2z`$ for the food webs and indeed find a collapse of the data onto a single curve for $`m/2z<0.7`$ (Figs. 2b,d) ks-test2 .

Figure 2 supports the strong new hypothesis that the distributions of number of prey and number of predators follow universal functional forms. To improve statistics and better determine the specific functional form of these distributions, one may pool the scaled variables, $`k/2z`$ and $`m/2z`$, from the different webs ks-test1 ; ks-test2  into single distributions, $`p_{\mathrm{prey}}`$ and $`p_{\mathrm{pred}}`$, respectively. Figures 3a,b show the cumulative distributions of scaled number of prey and scaled number of predators. Note that the distributions are well approximated by Eqs. (1)–(2) even though there are no free parameters to fit in the analytical curves. These results are analogous to the finding of scaling and universality in physical, chemical and social systems Bunde94 .

Figure 3c plots the probability densities for the distribution of number of prey and number of predators. It is visually apparent that both distributions are different. This is confirmed by the Kolmogorov-Smirnov test which rejects the null hypothesis at the $`p<0.001`$ level. The distribution of number of prey decays exponentially, and the distribution of number of predators is essentially a step function with a fast decay.

One can perform a similar analysis for the distribution $`p_{\mathrm{link}}`$ of the number of trophic links $`rk+m`$. As for number of prey or number of predators, the data from the different webs, upon the scaling $`r/2z`$, collapse onto a single curve, further supporting the hypothesis that scaling holds for food web structure. To better determine the specific functional form of $`p_{\mathrm{link}}(r)`$, we pool the scaled variables, $`r/2z`$, from all webs except Ythan into a single distribution (Fig. 3d). We find that $`p_{\mathrm{link}}(r)`$ has an exponential decay for $`r/2z1`$, in agreement with our theoretical calculations. Therefore, there is a characteristic scale for the linkage density, i.e. food webs do not have a scale-free structure, in contrast to reports in recent studies of food-web fragility Sole .

Next, we test if the scaling hypothesis suggested by the analysis of distribution of trophic links also applies to other quantities characterizing food web structure. We consider two quantities with ecologic implications: (i) the average trophic distance $`d`$ between species Watts98  (which is the number of species needed to trophically connect two given species), and (ii) the clustering coefficient $`C`$ (which counts the fraction of species’ triplets that form fully-connected triangles). The latter relates to the compartmentalization in an ecosystem while the former relates to typical food-chain length.

In Fig. 4a, we compare our numerical results for the average trophic distance $`d`$ for the niche model Martinez00  with the values calculated for the food webs analyzed. We find that $`d`$ increases with web size as $`\mathrm{log}S`$ both for the model and for the data. This logarithmic increase is the expected behavior for a random graph; however, the slopes measured for the data and the model are different from the value predicted for a random graph Newman01 , suggesting that there is a degree of “order” to the connectivity of the food web which may encode the mechanisms of food web assembly. Remarkably, this characteristic of the empirical food webs appears to be captured by the niche model Martinez00 . The results of Fig. 4a also support the scaling hypothesis and suggest that the average distance in a food web may also follow a unique functional form for different food webs.

Figure 4b shows our results for the clustering $`C`$ of the food webs studied and for the niche model Martinez00 . We find that the data is well approximated by the model predictions, and that $`C`$ decreases to zero as $`1/S`$ as web size $`S`$ increases.

The major finding of this study is the uncovering of unifying quantitative patterns characterizing the structure of food webs from diverse environments. Specifically, we find that the distributions of the number of prey, number of predators and number of links of most of the best studied food webs seem to collapse onto the same curves after rescaling the number of links by its average number $`z`$. Remarkably, the corresponding curves are in agreement with the analytical predictions of the niche model. Therefore, these distributions can be theoretically predicted merely by knowing the food web’s linkage density $`z`$, a parameter readily accessible empirically. Regularities such as these are interesting as descriptors of trophic interactions inside communities because they may enable us to make predictions in the absence of high-quality data, and provide insight into how communities function and are assembled.

Our results are of interest for a number of other reasons. First, food webs do not have a scale free distribution of number of links (total, incoming or outgoing). This is surprising since one could expect most species to try to prey on the most abundant species in the ecosystem (an “abundant-get-eaten” kind of mechanism). Such a preferential attachment would lead to a scale-free distribution of links; instead, we find a single-scale distribution, suggesting that species specialize and prey on a small set of other species. Second, the results of Figs. 4a,b support the scaling hypothesis and indicate that there is very little, if any, compartmentalization in ecosystems Pimm80 , suggesting the possibility that ecosystems are highly interconnected and that the removal of any species may induce large disturbances. Third, the structure of food webs is different from many other biological networks in two important aspects: the links are uni-directional and the in- and out-degree distribution are different. These two facts are a result of the directed character of the trophic interactions and of the asymmetry it creates. Interestingly, the niche model captures this asymmetry in its rules, which may explain its success in explaining the empirical results.

Our findings are surprising for two reasons: (i) they hold for the most complete food webs studied, in contrast to previously reported patterns ecology , and (ii) they support the possibility that fundamental concepts of modern statistical physics such as scaling and universality Bunde94 —which were developed for the study of inanimate systems—may also be applied in the study of food webs—which comprise animate beings. Indeed, our results are consistent with the underlying hypothesis of scaling theory Bunde94 , i.e., food webs display universal patterns in the way trophic relations are established despite apparently “fundamental” differences in factors such as the environment (e.g. marine versus terrestrial), ecosystem assembly, and past history.

We thank A. Arenas, J. Bafaluy, M. C. Barbosa, M. Barthélémy, A. Díaz-Guilera, J. Faraudo, G. Franzese, F. Giralt, S. Mossa, R. V. Solé, H. E. Stanley and especially N. D. Martinez and R. J. Williams for stimulating discussions and helpful suggestions. We also thank J. E. Cohen, N. D. Martinez and R. J. Williams for making their electronic databases of food webs available to us. JC and RG thank the Generalitat de Catalunya and the Spanish CICYT (PPQ2000-1339, BFM2000-0626 and BFM2000-0351-C03-01) for support. LANA thanks NIH/NCRR (P41 RR13622) and NSF for support.


text:
# Ab initio Hartree-Fock with electronic correlation study of the electronic properties of MgB2

## I INTRODUCTION

MgB2 has attracted lots of attention of the international scientific community due that it shows a superconducting transition of the order of 39 K. Furthermore, its unitary cell (three atoms per unit cell), is simple enough that provides a tool to understand the mechanism proposed for this material to behave like a BCS type of superconductor i.e. electron-phonon interaction.

Early theoretical studies of the electronic structure, using other methods of calculations like OPW and tight-binding, have shown similar behavior to graphite, from the point of view of band theory. Transition metal diborides (TMD), as expected, show a more complex structure due to the consideration of d-states.

Recently, an electronic structure calculation on MgB2 shows that some states of boron are responsible for the metallic behavior of this compound. Furthermore, phonon frequencies and electron-phonon coupling (at the most symmetric point of the Brillouin zone) lead these authors to conclude to phonon based superconductivity.

## II  COMPUTATIONAL DETAILS

The calculations reported in the present work had been carried out using CRYSTAL98 ab initio Self-consistent-field (SCF) Hartree-Fock linear combination of atomic orbital computer program which provide solution to crystalline system of any space group. Electronic correlation using density functional (DFT) has been included as a correction to the total energy using the correlation function proposed by Perdew-Zunger (PZ). More details about the mathematical formulation of CRYSTAL98 have been described else where and will be omitted here.

The crystal structure of MgB2 is illustrated in Fig. 1. It has a hexagonal unit cell, with the primitive vectors $`a=3.086`$ Å and $`c=3.504`$ Å respectively. The space group for this compound is P6/mmm (No. 191). From the figure is possible to notice that the boron atoms are arranged in layers, with the next layers of Mg interleaved between them. In order to form the unit cell two layers of magnesium and one of boron are needed along the c-axis in the hexagonal configuration. For the basis set for magnesium and boron, we selected 6-21G\* and 6-21G basis set respectively, provided in Ref. (11). The selection of the basis set was obtained between a compromise between accuracy and computational resources.

## III  RESULTS AND DISCUSSION

For this calculation, we used the experimental lattice parameters provided by Nagamatsu et al. of a= 3.086 Å, c= 3.524 Å respectively. Band structure calculations for MgB2 are depicted in Fig. 2. In order to obtain the band structure for this compound, we used 133 k-points sampling the FBZ.

Our results are in good qualitative agreement with the results obtained by other theoretical methods, such as the results obtained by Kortus et al. who used a general potential LAPW code, whilst Satta et al.  using local-density approximation to the density functional theory, and full-potential linearized augmented plane waves, obtained similar results. Furthermore, Belaschenko calculated the electronic structure for this compound using Stuttgard LMTO-TB (ASA) to compute the energy bands with similar results.

We spanned the FBZ covering the three dimensional space going from $`\mathrm{\Gamma }`$-M-K-$`\mathrm{\Gamma }`$-A-L-H-A of $`k_x=p/a,k_y=p/b,k_z=p/c`$. The $`\mathrm{\Gamma }`$-M-K-$`\mathrm{\Gamma }`$ lines are in the basal plane, while the A-L-H-A are located on the top of the plane at $`k_z`$. The upper part of the valence band for MgB2, composed of B 2p-states which form two peculiar set of bands of $`\sigma `$(2px,y) and $`\pi `$(pz) character, whose $`k`$ dependence changes considerably while spanning the three dimensional space. For B 2px,y the most pronounced dispersion is along $`\mathrm{\Gamma }`$-K. These bands are of quasi-two dimensional type and form an almost flat zone along $`\mathrm{\Gamma }`$-A, which reflects the distribution of pp$`\sigma `$-states in the B layers. These states provide a big contribution to the density of states close to the Fermi energy, providing the metallic properties to the borides.

On the other hand, B pz like bands are responsible to the weak pp$`\pi `$ interactions. These bands are considered to be a three dimensional like type, which have maximum dispersion along $`\mathrm{\Gamma }`$-A. Furthermore, Mg s, p and B s states are admixed with B 2p bands close to the bottom of the valence band, as well as in the conduction band. Hence the electronic properties of MgB2 are associated with metallic 2p-states of B atoms located in planar nets, which determines the DOS close to the Fermi level.

For the total and projected density of states (PDOS) for B and Mg, see Fig. 3. Concentrating on those states close to the Fermi level, which presumably are responsible for the superconductor behavior of the material, as we mentioned before in the energy band discussion, the main contributions comes from p-$`\sigma `$ bonding and pz having $`\pi `$ bonding and anti bonding character on the basal and top planes, respectively. Which one of these two bands is responsible for the superconductor behavior is unclear so far.

In summary, we have shown energy bands, as well as total and projected density of states for the medium temperature superconductor MgB2. From these analysis, we were able to infer the most important contributions from each atom to the superconductor behavior for the compound. So far, this study is not conclusive, but continues in order to provide other important physical properties such as charge density profiles and Mulliken population analysis.

###### Acknowledgements.

Both of us acknowledge DGSCA-UNAM and CoNaCyT for providing financial support.


text:
# Comment on studying the corrections to factorization in 𝐵→𝐷^(∗)𝑋

## Abstract

We propose studying the mechanism of factorization in exclusive decays of the form $`BD^{()}X`$ by examining the differential decay rate as a function of the invariant mass of the light hadronic state $`X`$. If factorization works primarily due to the large $`N_c`$ limit then its accuracy is not expected to decrease as the $`X`$ invariant mass increases. However, if factorization is mostly a consequence of perturbative QCD then the corrections should grow with the $`X`$ invariant mass. Combining data for hadronic tau decays and semileptonic $`B`$ decays allows tests of factorization to be made for a variety of final states. We discuss the examples of $`BD^{}\pi ^+\pi ^{}\pi ^{}\pi ^0`$ and $`BD^{}\omega \pi ^{}`$. The mode $`BD^{}\omega \pi ^{}`$ will allow a precision study of the dependence of the corrections to factorization on the invariant mass of the light hadronic state.

preprint:   LBNL–47551 UTPT–01–05 CALT–68–2316 hep-ph/0103020

An understanding of nonleptonic $`B`$ decays is important for the study of $`CP`$ violation at $`B`$ factories. However, there is at present no systematic, model-independent treatment of all such decays. In certain kinematic situations, factorization, according to which the matrix element of four-quark operators can be written as the product of pairs of matrix elements of two-quark operators, can be justified using perturbative QCD . An example is $`BD^{()}X`$, where $`X`$ denotes a hadronic state with low invariant mass, and the $`D^{()}`$ inherits the light degrees of freedom from the $`B`$ meson. The large $`N_c`$ limit of QCD  provides a different explanation of factorization, which is independent of the kinematics of the final state. Factorization has been shown to be consistent with experiment in certain two-body decays, such as $`BD^{()}\pi `$ and $`D^{()}\rho `$. Currently, this represents a test of factorization at the $`10`$$`20\%`$ level . Factorization has also been observed to hold within the (presently sizable) errors in decays of the type $`BD^{()}D_s^{()}`$, where no perturbative QCD justification has been presented .

In this paper we examine how the decays $`BD^{()}X`$, where $`X`$ is a final state containing two or more light hadrons, allow us to study the corrections to factorization***These decays have also been discussed as tests of factorization by Reader and Isgur., and help determine the role that perturbative QCD plays. In such decays, factorization can be studied as a function of $`m_X`$, the invariant mass of the final hadronic state produced by the light quark current. The perturbative QCD arguments for factorization depend on the light quarks being produced in an almost collinear state (for this reason, it has also been argued that factorization should hold for $`BD^{()}+\mathrm{jet}`$). Corrections are therefore expected to grow with $`m_X/E_X`$, where $`E_X`$ is the energy of the hadronic final state in the $`B`$ rest frame, since this quantity characterizes the deviation of the hadronic state $`X`$ from the light cone. However, if factorization works primarily due to the large $`N_c`$ limit then its accuracy is not expected to decrease as $`m_X`$ increases. The behavior of factorization violations with $`m_X`$ can therefore distinguish the perturbative QCD explanation of factorization from that of the large $`N_c`$ limit. Naïvely, the $`1/N_c^2`$ terms that violate factorization are $`O(10\%)`$ in the amplitude and $`O(20\%)`$ in the rate. When factorization is tested with a greater precision than this, it may be possible to see corrections to factorization which grow with $`m_X`$, which would provide evidence that perturbative QCD plays an important role in factorization.

Nonleptonic $`B`$ meson decays arise dominantly due to the weak Hamiltonian

$$H^{\mathrm{nl}}=\frac{4G_F}{\sqrt{2}}V_{cb}V_{uq}^{}\left[c_1(m_b)(\overline{c}\gamma _\mu P_Lb)(\overline{q}\gamma ^\mu P_Lu)+c_2(m_b)(\overline{q}\gamma _\mu P_Lb)(\overline{c}\gamma ^\mu P_Lu)\right],$$
(1)
where $`P_L=(1\gamma _5)/2`$ and $`q=d,s`$. The Hamiltonian in Eq. (1) is renormalized at the scale $`\mu =m_b`$. The coefficients $`c_i`$ are given by $`c_{1,2}=(c_+\pm c_{})/2`$, where

$$c_+(\mu )=\left[\frac{\alpha _s(\mu )}{\alpha _s(M_W)}\right]^{6/23},c_{}(\mu )=\left[\frac{\alpha _s(\mu )}{\alpha _s(M_W)}\right]^{12/23},$$
(2)
in the leading logarithmic approximation. At next-to-leading order, $`c_1(m_b)=1.13`$ and $`c_2(m_b)=0.29`$ for $`\alpha _s(M_Z)=0.118`$. Under the factorization hypothesis, the hadronic matrix element of $`H^{\mathrm{nl}}`$ can be written as

$$XD^{()}|H^{\mathrm{nl}}|B=\frac{4G_F}{\sqrt{2}}V_{cb}V_{uq}^{}\left(c_1(m_b)+\frac{1}{3}c_2(m_b)\right)D^{()}|\overline{c}\gamma _\mu P_Lb|BX|\overline{q}\gamma ^\mu P_Lu|0+\mathrm{},$$
(3)
where the ellipses denote factorization violating terms. The first matrix element can be determined from semileptonic $`BD^{()}\mathrm{}\overline{\nu }`$ decays, while the second can be determined from the hadronic tau decay $`\tau X\nu _\tau `$.

The differential rate for tau decay to a given hadronic final state $`X`$ may be written as

$$\frac{\mathrm{d}\mathrm{\Gamma }(\tau X\nu _\tau )}{\mathrm{d}m_X^2}=\frac{G_F^2|V_{ud}|^2}{32\pi ^2m_\tau ^3}(m_\tau ^2m_X^2)^2(m_\tau ^2+2m_X^2)v_X(m_X^2),$$
(4)
where $`m_X^2`$ is the invariant mass-squared of the hadronic final state $`X`$, and $`v_X(m_X^2)`$ is related to the transverse part of the vacuum polarization amplitude. This holds in general even if $`X`$ is a nonresonant multibody final state. The only assumption is that the weak hadronic current is conserved. For the vector part of the current, corrections to this assumption are suppressed by $`m_um_d`$, while for the axial part they are suppressed by $`m_u+m_d`$. Perturbative electroweak corrections have been neglected in Eq. (4).

The CLEO collaboration has measured $`v_X(m_X^2)`$ for two-, three-, and four-pion final states. The two-pion final states are quasi-two-body, being dominated by the $`\rho `$ resonance. Since the $`\rho `$ is narrow, $`\mathrm{d}\mathrm{\Gamma }(BD^{()}\pi \pi )/\mathrm{d}m_X^2`$ is a steeply falling function of $`m_X^2`$ in the region we are interested in, and so is less useful for our purposes than the three- and four-pion final states. The three-pion final states receive a large $`a_1`$ contribution, however, the $`a_1`$ is rather broad and there is evidence for nonresonant contributions in $`\tau `$ decay to three pions. The four-pion final states are not dominated by any one resonance.

In terms of $`v_X(m_X^2)`$, the factorization prediction for the decay $`BD^{()}X`$ is

$$\frac{\mathrm{d}\mathrm{\Gamma }(BD^{()}X)/\mathrm{d}m_X^2}{\mathrm{d}\mathrm{\Gamma }(BD^{()}\mathrm{}\overline{\nu })/\mathrm{d}m_X^2}=3\pi \left(c_1(m_b)+\frac{c_2(m_b)}{3}\right)^2v_X(m_X^2)(1+\delta _{NF}),$$
(5)
where $`\delta _{\mathrm{NF}}`$ denotes the nonfactorizable contributions, and in the denominator of the left-hand side $`m_X`$ denotes the dilepton invariant mass. In the large $`N_c`$ limit $`\delta _{NF}`$ is of order $`1/N_c`$, whereas in the perturbative QCD approach it is expected to contain both perturbative corrections and nonperturbative ones suppressed by powers of $`m_X/E_X`$, where $`E_X=(m_B^2+m_X^2m_{D^{()}}^2)/(2m_B)`$ is the energy of the state $`X`$ in the $`B`$ rest frame. The ratio $`m_X/E_X`$ changes from $`0.43`$ at $`m_X=1`$GeV to $`0.70`$ at $`m_X=m_\tau `$, so if there were order $`m_X/E_X`$ corrections to factorization then the accuracy of the predictions should change significantly over the accessible range of $`m_X`$.

In Fig. 1 we plot the $`m_X^2`$ distribution in $`BD^{}\pi ^+\pi ^{}\pi ^{}\pi ^0`$ decay, normalized to the $`BD^{}\mathrm{}\overline{\nu }`$ rate, as predicted by Eq. (5). We used the CLEO fit to the shape of $`\mathrm{d}\mathrm{\Gamma }(BD^{}\mathrm{}\overline{\nu })/\mathrm{d}m_X^2`$  and the $`\tau `$ decay data . This decay mode is interesting for testing factorization because it is not dominated by a single narrow resonance, and its branching fraction is large, almost $`2\%`$ . However, this mode has a significant disadvantage for testing factorization because of a potentially sizeable background from processes where one or more of the pions are emitted from the $`(\overline{c}b)`$ current. In this case, the $`(\overline{c}b)`$ current creates a charm state $`C`$ of the type $`D^{}+n\pi `$ ($`1n3`$), and the $`(\overline{d}u)`$ current creates the remaining $`4n`$ pions in the final state. In the formal limit when $`m_bm_c`$ is large this background is parametrically suppressed, since the light hadrons from the $`C`$ decay combined with those from the $`(\overline{d}u)`$ current will have large invariant mass over most of the phase space. This can, however, be a sizable background for the physical $`m_b`$ and $`m_c`$ masses.

Because of the problem of backgrounds of this type, a more promising way to test the $`m_X^2`$ dependence of factorization is to focus on the subset of these decays which proceeds via $`BD^{}\omega \pi ^{}`$. In this case, the backgrounds are easier to constrain, because the $`(\overline{d}u)`$ current must create the $`\pi ^{}`$ in the final state. Therefore, the $`(\overline{c}b)`$ current would have to create charmed states, $`C_{(D^{}\omega )}`$, which can decay into $`D^{}\omega `$. This background can be related using factorization (in the most reliable case of the current only making one pion) to the $`BC_{(D^{}\omega )}\mathrm{}\overline{\nu }`$ decay rate. Assuming that $`C_{(D^{}\omega )}`$ is a narrow resonance,

$$\mathrm{\Gamma }(BC\pi )=\frac{3\pi ^2|V_{ud}|^2f_\pi ^2}{m_Bm_C}\left(c_1(m_b)+\frac{c_2(m_b)}{3}\right)^2\left(\frac{\mathrm{d}\mathrm{\Gamma }(BC\mathrm{}\overline{\nu })}{\mathrm{d}w}\right)_{w_{\mathrm{max}}},$$
(6)
where $`w=v_Bv_C`$. Assuming that the $`C_{(D^{}\omega )}`$ mass is its minimal value $`m_C^{\mathrm{min}}=m_D^{}+m_\omega `$ and that semileptonic $`BC_{(D^{}\omega )}\mathrm{}\overline{\nu }`$ rate depends on $`w`$ as $`\sqrt{w^21}`$, we find that $`\mathrm{\Gamma }(BC_{(D^{}\omega )}\pi )0.1\times \mathrm{\Gamma }(BC_{(D^{}\omega )}\mathrm{}\overline{\nu })`$. Therefore, this background can be constrained by measuring the $`BD^{}\omega \mathrm{}\overline{\nu }`$ decay rate. The measured $`BD^{}\omega \pi ^{}`$ branching fraction is about $`0.3\%`$ . Thus, for example, a bound of $`(BC_{(D^{}\omega )}\mathrm{}\overline{\nu })<0.3\%`$ would imply that the feed-down from higher mass $`C_{(D^{}\omega )}`$ states to the $`D^{}\omega \pi ^{}`$ final state is less than $`10\%`$. This is a very conservative estimate, since the $`C_{(D^{}\omega )}`$ branching fraction to $`D^{}\omega `$ is probably significantly less than unity. Moreover, the $`BC_{(D^{}\omega )}\mathrm{}\overline{\nu }`$ branching fraction is probably much smaller than $`0.3\%`$; this is certainly the case in most models. Therefore, we expect that the $`BD^{}\omega \pi ^{}`$ decay will only receive a very small background from higher mass charm states.

In Fig. 2 we plot the $`m_X^2`$ distribution in $`BD^{}\omega \pi ^{}`$ decay, normalized to the $`BD^{}\mathrm{}\overline{\nu }`$ rate, as predicted by Eq. (5) and the $`\tau `$ decay data . This decay receives a large contribution from a $`\rho ^{}`$ resonance with mass around $`1.41.5`$GeV . However, the $`\rho ^{}`$ is very broad, with a width estimated to be around $`400500`$MeV. Since we would like to make a precision test of factorization and learn about the size of corrections suppressed by powers of $`m_X/E_X`$ and $`\alpha _s`$, one needs tests that do not make assumptions based on resonance models or heavy quark symmetry , which may have corrections of similar size. This is best achieved by looking at the differential mass distribution, which should also teach us about how factorization breaks down (if it does at all) at higher masses.

Controlling this background in the $`BD^{}\pi ^+\pi ^{}\pi ^{}\pi ^0`$ decay is more complicated than for $`BD^{}\omega \pi ^{}`$, and may make it difficult to make a precise study of factorization in this channel. The contribution when $`C`$ contains two or three pions is expected to be small and can be constrained by bounding the semileptonic rates as in the case of $`BD^{}\omega \pi ^{}`$. The contribution when one pion comes from the $`(\overline{c}b)`$ current and three pions from the $`(\overline{d}u)`$ current may be estimated from the measured $`v_{3\pi }(m_X^2)`$ and the $`BD^{}\pi \mathrm{}\overline{\nu }`$ branching fraction. Using the ALEPH data  for $`v_{3\pi }(m_X^2)`$ measured in $`\tau `$ decay and $`(BD^{}\pi \mathrm{}\overline{\nu })1\%`$, it appears that this background may be significant. The size of this background can be tested experimentally by measuring the $`BD^0\pi ^+\pi ^+\pi ^{}\pi ^{}`$ decay rate, since this decay only receives contributions from higher mass charm states which decay into a $`D^0`$ plus one or three pions. Unless the $`BD^0\pi ^+\pi ^+\pi ^{}\pi ^{}`$ rate is much smaller than the $`BD^{}\pi ^+\pi ^{}\pi ^{}\pi ^0`$ rate, it will be difficult to make a precise test of factorization in $`BD^{}\pi ^+\pi ^{}\pi ^{}\pi ^0`$.

Data on $`e^+e^{}`$hadrons can also be used to determine $`v_X(m_X^2)`$ and predict the mass spectrum in $`BD^{()}X`$ decay . Currently the CMD–2 data on $`e^+e^{}4\pi `$ are available up to $`m_X^2=(1.4\mathrm{GeV})^2`$ ; measurements of this cross section at higher energies would allow testing factorization at values of $`m_X^2>m_\tau ^2`$.

In summary, for a precision study of corrections to factorization in $`BD^{()}X`$ decays (where $`X`$ contains only light hadrons), it is best not to make assumptions based on heavy quark symmetry or resonant decomposition (especially for broad states). Rather, one should use the differential $`m_X^2`$ distribution with a fixed hadron content and compare with the predicted rate based on factorization using the spectral function extracted from $`\tau `$ decay. Over the accessible kinematic range the parameter $`m_X/E_X`$, which controls the accuracy of the predictions if the primary reason for factorization is perturbative QCD, changes from about $`0.4`$ to $`0.7`$, making $`m_X`$-dependent violations of factorization potentially observable. In this paper we concentrated on final states containing four pions, but this approach is general and can be applied to test factorization in decays to any light hadronic final state where backgrounds from the $`(\overline{c}b)`$ current producing some of the light hadrons can be controlled. Studying the available data in $`BD^{}\pi ^+\pi ^{}\pi ^{}\pi ^0`$ and $`BD^{}\omega \pi ^{}`$ decays, we found no evidence of factorization becoming a worse approximation at higher values of $`m_X`$. It will be very interesting to make a precision study, which should be possible at the $`B`$ factories. The decay mode $`BD^{}\omega \pi ^{}`$ is particularly well suited for this purpose.

###### Acknowledgements.

We thank Sheldon Stone, Jianchun Wang and Alan Weinstein for help with the CLEO data and discussions. Z.L. was supported in part by the Director, Office of Science, Office of High Energy and Nuclear Physics, Division of High Energy Physics, of the U.S. Department of Energy under Contract DE-AC03-76SF00098. M.L. was supported in part by the Natural Sciences and Engineering Research Council of Canada. M.B.W. was supported in part by the Department of Energy under Grant No. DE-FG03-92-ER40701.