text:
# Improved Measurement of |𝑉_{𝑢𝑏}| with Inclusive Semileptonic 𝐵 Decays

CLEO Collaboration

(February 7, 2002)

## Abstract

We report a new measurement of the CKM parameter $`|V_{ub}|`$ made with a sample of 9.7 million $`B\overline{B}`$ events collected with the CLEO II detector. Using Heavy Quark theory, we combine the observed yield of leptons from semileptonic $`B`$ decay in the end-point momentum interval $`2.22.6`$ GeV/$`c`$ with recent CLEO II data on $`BX_s\gamma `$ to find $`|V_{ub}|=(4.08\pm 0.34\pm 0.44\pm 0.16\pm 0.24)\times 10^3`$, where the first two uncertainties are experimental and the last two are from theory.

preprint:   CLNS 01/1767 CLEO 01-23

Measurement of the Cabibbo-Kobayashi-Maskawa matrix element $`V_{ub}`$ is one of the most challenging and important tasks in testing the Standard Model. This direct measurement of one side of the unitarity triangle gives information that is complementary to studies of CP violation in $`B`$ decays. Semileptonic $`B`$ decays are a direct route to $`|V_{ub}|`$ and $`|V_{cb}|`$. Measurements of lepton production above the kinematic limit for $`BX_c\mathrm{}\nu `$ provided the discovery of nonzero $`|V_{ub}|`$ Fulton et al. (1990); Albrecht et al. (1990) and a measurement of its value Bartelt et al. (1993) that was controversial due to its heavy reliance on theoretical models.

In this paper we report an improved inclusive measurement of $`BX_u\mathrm{}\nu `$ made with the CLEO II detector at the Cornell Electron Storage Ring (CESR). This measurement supersedes that of Ref. Bartelt et al. (1993), which was based on the first tenth of our data sample. Improved background-suppression techniques help to reduce the model dependence of the result. We have also measured $`BX_u\mathrm{}\nu `$ over a broader momentum range than was previously possible, further helping to reduce theoretical uncertainty, but incurring experimental uncertainty because of the large $`BX_c\mathrm{}\nu `$ subtraction. We use the momentum range $`2.22.6`$ GeV/$`c`$ as a compromise between these considerations, and present other intervals to probe the stability of our result and to compare with past measurements. For extracting $`|V_{ub}|`$, models have been replaced with Heavy Quark (HQ) theory Neubert (1994a, b); Bigi et al. (1994a, b); Kagan and Neubert (1999); De Fazio and Neubert (1999); Leibovich et al. (2000) and the CLEO-measured $`BX_s\gamma `$ photon-energy spectrum Chen et al. (2001). These improvements result in a determination of $`|V_{ub}|`$ with overall uncertainty that is smaller and more reliably estimated than previous measurements.

The CLEO II detector has been described in detail Kubota et al. (1992). Key components for this work are the tracking system, the cesium iodide calorimeter, and the muon detector. Two thirds of our data were collected after a detector upgrade that included a silicon tracker Hill (1998), a change of drift-chamber gas from argon-ethane to helium-propane, and other improvements. We obtained an integrated luminosity of 9.13 fb-1 of $`e^+e^{}`$ annihilation data at the $`\mathrm{{\rm Y}}(4S)`$ resonance (ON) and 4.35 fb-1 at energies just below $`B\overline{B}`$ threshold (OFF). OFF yields are scaled by a factor of $`2.06\pm 0.02`$ and subtracted from the ON yields to eliminate contributions of non-$`B\overline{B}`$ processes. We determine the scale factor by comparing ON and OFF track spectra above the $`B`$-decay kinematic limit, and confirm the result with measured luminosities. The ON sample includes 9.7 million $`B\overline{B}`$ events.

We select electron and muon candidates in the central region of CLEO II ($`|\mathrm{cos}\theta |0.71`$, where $`\theta `$ is the angle with respect to the beam). An electron candidate must have an energy deposit in the calorimeter close in value to its measured momentum and a specific ionization consistent with that expected for an electron. The efficiency of electron selection is slightly momentum-dependent, with a value of $`0.93\pm 0.03`$ at 2.2 GeV/$`c`$. Muons must penetrate seven nuclear interaction lengths of absorber material. The threshold momentum is $`1.8`$ GeV/$`c`$ and the efficiency is $`0.88\pm 0.03`$ at 2.2 GeV/$`c`$. Misidentification rates are measured for tagged $`\pi ^\pm `$, $`K^\pm `$, and $`p/\overline{p}`$ tracks and averaged in proportions determined with simulated $`B\overline{B}`$ events. The electron and muon fake rates above 2.0 GeV/$`c`$ are less than 0.1% and about 0.5%, respectively.

We use several cuts to minimize mismeasurements and background from sources other than semileptonic decays. We veto any lepton that can be paired with another lepton of the same type and opposite charge if the pair mass is within three standard deviations of the $`J/\psi `$ mass. Electrons from photon conversions are also rejected. Stringent requirements on tracking residuals, impact parameter, and the fraction of tracking layers traversed that have good hits ensure a reliable momentum measurement. The momentum-dependent track-selection efficiency is determined with a GEANT simulation Brun et al. (1987) of signal events. Real and simulated radiative Bhabha electrons embedded in hadronic events are used to determine a correction factor accounting for the imperfect simulation. The track-selection efficiency for $`BX_u\mathrm{}\nu `$ is $`0.95`$, with a 3% uncertainty.

Although most leptons at the $`\mathrm{{\rm Y}}(4S)`$ are from semileptonic $`B`$ decays, continuum production contributes a large background at high momentum that can be subtracted with OFF data. To avoid degrading the precision of our $`BX_u\mathrm{}\nu `$ measurement, we use several criteria to suppress the continuum before subtraction. Beam-gas, beam-wall, and some QED events are eliminated by cuts on visible energy and the event vertex. A multiplicity cut suppresses $`\tau ^+\tau ^{}`$. Particles above the $`B`$-decay kinematic limit indicate continuum production or mismeasurement, so we eliminate events having a track or shower with an energy above 3.5 GeV. Two-photon and QED processes are suppressed by cutting events with missing momentum near the beam.

The continuum background remaining after this selection is primarily $`e^+e^{}c\overline{c}`$. In our previous analysis Bartelt et al. (1993), cuts on event-shape variables and missing momentum used the kinematics of $`BX_u\mathrm{}\nu `$ to suppress the continuum by a factor of 70, with a signal efficiency of $`0.38`$. The efficiency depended strongly on $`q^2`$ (the squared mass of the virtual $`W`$), and thus contributed to the model dependence of the result. The goal of this analysis is to achieve comparable suppression with reduced model dependence. This requires cuts that exploit the presence of the “other $`B`$” in signal events rather than the details of $`BX_u\mathrm{}\nu `$. We use a neural net with inputs that are measurements of the energy flowing into eleven angular intervals defined with respect to the candidate-lepton direction. We exclude energy within 25 of the direction opposite to the lepton, where there is strong dependence on the $`q^2`$ of the decay. The net was trained with a simulated signal sample generated with the ISGW2 model Scora and Isgur (1995) and a simulated background sample of continuum $`e^+e^{}`$ $`q\overline{q}`$. Optimized for expected signal level, the neural net gives background rejection of a factor of 50 and signal efficiency of 0.33.

The total efficiency for selecting leptons from $`BX_u\mathrm{}\nu `$ ($`0.21`$) shows a small momentum dependence over the range $`2.02.6`$ GeV/$`c`$. The uncertainty is about 7%, with roughly equal contributions from detector response and model uncertainties. We assess the latter by considering several models: ISGW2, the ACCMM spectator model Altarelli et al. (1982); Artuso (1993) with various values of internal parameters, and a hybrid of these developed by CLEO. Our current procedures have less than one third of the model-to-model variation of the previous analysis and essentially no dependence on the mass of the hadronic system $`X_u`$.

The computation of the $`BX_u\mathrm{}\nu `$ signal for the momentum interval $`2.22.6`$ GeV/$`c`$ is shown in Table 1. The spectrum of lepton candidates from $`B\overline{B}`$ events is obtained by subtracting the scaled OFF momentum distribution from the ON distribution. It includes $`BX_u\mathrm{}\nu `$ and several $`B`$-decay backgrounds. Hadrons misidentified as leptons (fakes) are computed by combining the momentum-dependent misidentification probabilities with the spectrum of hadrons from $`B\overline{B}`$ that pass all selection requirements except for lepton identification. Leakage through the $`J/\psi `$ and photon-conversion vetoes is estimated with Monte Carlo normalized to data, with systematic errors of 10% and 25%, respectively. Several other backgrounds are also estimated by Monte Carlo and are assigned conservative systematic errors: $`\psi (2S)`$ leptonic decays ($`\pm 25\%`$), semileptonic and leptonic decays of $`D`$ and $`D_s`$ ($`\pm 50\%`$), semileptonic decays to $`\tau `$ ($`\pm 25\%`$), and $`\pi ^0`$ Dalitz decays ($`\pm 100\%`$).

The largest background to $`BX_u\mathrm{}\nu `$ for momenta below 2.4 GeV/$`c`$ is $`BX_c\mathrm{}\nu `$. We calculate this by fitting the lepton spectra between 1.5 GeV/$`c`$ and the low end of the end-point interval. Fit functions are generated with Monte Carlo simulations using models and CLEO-measured form factors Duboscq et al. (1996); Athanas et al. (1997). QED radiative corrections are applied with the PHOTOS algorithm Barberio and Wa̧s (1994). Spectra generated in the $`B`$ rest frame are boosted to the lab frame using the $`B`$-momentum distribution of our data. There are four components in the fits: a mixture of $`D`$ and $`D^{}`$ with the ratio given by exclusive branching fractions Groom et al. (2000); a mixture of decays to $`D^{}`$ and other higher-mass charmed mesons (ISGW2); nonresonant decays (model of Goity and Roberts Goity and Roberts (1995)); and $`BX_u\mathrm{}\nu `$ (ISGW2 normalized to the end-point yield of Ref. Bartelt et al. (1993)).

The proportions for the $`BX_c\mathrm{}\nu `$ components are determined by fits to the electron spectrum without the neural-net cut, because the momentum acceptance is larger and has less uncertainty than that for the muons. The resulting fit has a $`\chi ^2`$ of 14.9 for 11 degrees of freedom. The muons are then fitted to the same mixture, with one parameter to allow for a difference in the $`e/\mu `$ efficiency ratio between Monte Carlo and data. The electron and muon spectra with the neural-net cut applied are then fitted to the same mixtures, with one parameter allowing for mis-modeling of the neural-net efficiency. These fits provide the $`BX_c\mathrm{}\nu `$ subtractions in Table 1.

We assess the uncertainty in the $`BX_c\mathrm{}\nu `$ subtraction by varying inputs, including the $`D/D^{}`$ ratio, the form factors for $`BD/D^{}\mathrm{}\nu `$, the relative $`D^{}`$ and nonresonant normalizations, the radiative corrections, the $`B`$-momentum scale, the normalization and models for $`BX_u\mathrm{}\nu `$, and the fit intervals. The overall uncertainty in the subtraction is the sum in quadrature of the observed variations. The largest is due to form factors.

Fig. 1(a) shows the ON and scaled OFF momentum spectra.

Fig. 1(b) shows the background-subtracted and efficiency-corrected spectrum for $`BX_u\mathrm{}\nu `$. Statistical and systematic uncertainties have been combined in quadrature. Below 2.3 GeV/$`c`$, the $`BX_c\mathrm{}\nu `$ subtraction dominates the uncertainties, which are strongly correlated from bin to bin. The partial branching fraction for $`BX_u\mathrm{}\nu `$ ($`\mathrm{}`$ = $`e`$ or $`\mu `$) in a given momentum interval is $`\mathrm{\Delta }_u(p)=N_{\mathrm{}}(p)/[2N_{B\overline{B}}ϵ(p)]`$, where $`N_{\mathrm{}}`$ is half of the total yield of $`e`$’s and $`\mu `$’s, $`ϵ(p)`$ is the efficiency, and $`N_{B\overline{B}}=(9.67\pm 0.17)\times 10^6`$ is the number of $`B\overline{B}`$ events. Results are given in Table 2 for five different end-point momentum intervals.

They are in good agreement with our previous measurement Bartelt et al. (1993).

To determine the charmless semileptonic branching fraction $`(BX_u\mathrm{}\nu )`$ from the partial branching fraction $`\mathrm{\Delta }_u(p)`$, we need to know the true fraction $`f_u(p)`$ of the $`BX_u\mathrm{}\nu `$ spectrum that falls in the given momentum interval. Parton-level decays can be reliably calculated, but the observable meson-decay processes depend on the mass of the $`b`$ quark and its motion inside the $`B`$ meson. These properties have been described in recent years with HQ theory Neubert (1994a, b); Bigi et al. (1994a, b); Kagan and Neubert (1999); De Fazio and Neubert (1999); Leibovich et al. (2000). A light-cone shape function can be convoluted with the parton-level $`bu\mathrm{}\nu `$ to obtain the spectrum for $`BX_u\mathrm{}\nu `$. The shape function depends on nonperturbative QCD and has not yet been calculated from first principles. To leading order, the same shape function describes all $`b`$-to-light transitions, in particular also relating $`BX_s\gamma `$ to $`bs\gamma `$.

We have recently measured the photon-energy spectrum in $`BX_s\gamma `$ Chen et al. (2001). Three two-parameter descriptions of the shape function Bigi et al. (1994b); Kagan and Neubert (1999) are used to fit this spectrum over the range $`1.5<E_\gamma <2.8`$ GeV. All shape functions give good fits. We use the best-fit parameters and error ellipses to calculate the $`BX_u\mathrm{}\nu `$ lepton-momentum spectra and determine $`f_u`$ and its uncertainty, following Ref. De Fazio and Neubert (1999). We compute the effect of QED radiative corrections on $`f_u`$ with PHOTOS. The resulting values of $`f_u`$ are given in Table 2. The systematic uncertainty includes contributions from the subtraction of $`B`$-decay processes other than $`BX_s\gamma `$, the choice of scale for evaluating $`\alpha _s`$, radiative corrections, and differences among the shape functions.

We extract $`|V_{ub}|`$ by averaging the nearly identical formulations of Hoang et al. Hoang et al. (1999) and Uraltsev Uraltsev (1999):

$$|V_{ub}|=(3.07\pm 0.12)\times 10^3\times \left[\frac{(BX_ue\nu )}{0.001}\frac{1.6\mathrm{ps}}{\tau _B}\right]^{\frac{1}{2}}.$$

For the $`B`$ lifetime we use $`\tau _B=1.60\pm 0.02`$ ps Groom et al. (2000). Results for $`|V_{ub}|`$ (Table 2) show excellent agreement among the five momentum ranges. The best overall precision (15%) is obtained for the $`2.22.6`$ GeV/$`c`$ interval. The first uncertainty is from the measurement of $`\mathrm{\Delta }_u(p)`$ (combined statistical and systematic). The second is the combined uncertainty in the determination of $`f_u`$ from $`BX_s\gamma `$. The third is the average of the uncertainties given in Refs. Hoang et al. (1999); Uraltsev (1999) for the extraction of $`|V_{ub}|`$ from the branching fraction. The fourth is the theoretical uncertainty associated with the assumption that $`BX_s\gamma `$ can be used to compute the spectrum for $`BX_u\mathrm{}\nu `$. This is valid to leading order, with corrections at order $`\mathrm{\Lambda }_{QCD}/M_B`$. Taking $`\mathrm{\Lambda }_{QCD}/M_B0.1`$, we estimate this error by varying the shape function parameters by $`\pm `$10% Neubert . Uncertainty due to the theoretical assumption of quark-hadron duality remains unquantifiable.

For comparison, we also determine $`f_u`$ with models of $`BX_u\mathrm{}\nu `$. Available exclusive models are limited in the final states included, while inclusive models have uncertain internal parameters. Neither type provides a solid basis for assessing the theoretical uncertainty. For ISGW2 and ACCMM (spectator mass $`m_{sp}`$=150 MeV/$`c^2`$, Fermi momentum $`p_F`$=300 MeV/$`c`$), we find $`|V_{ub}|`$ in the $`2.22.6`$ GeV/$`c`$ interval to be 20% smaller than the result given above. Replacing the default ACCMM parameters (chosen for consistency with past analyses) with values obtained by fitting the $`BX_s\gamma `$ spectrum to a spectator parameterization Ali and Greub (1995) gives better agreement. The parameter values are $`m_{sp}`$230 MeV/$`c^2`$ and $`p_F`$440 MeV/$`c`$, leading to a $`|V_{ub}|`$ that is $``$5% smaller than our result.

In conclusion, we have measured the CKM parameter $`|V_{ub}|`$ to be $`(4.08\pm 0.63)\times 10^3`$. This result has smaller overall uncertainty than previous measurements and represents a major step forward, both in the quality of the experimental data and in the use of QCD theory rather than phenomenological models.

We gratefully acknowledge the efforts of the CESR staff in providing our excellent data sample. We thank A. Falk, A. Kagan, A. Leibovich, M. Luke, M. Neubert, I. Rothstein, M. Shifman, A. Vainshtein, and M. Voloshin for useful correspondences and discussions. This work was supported by the National Science Foundation, the U.S. Department of Energy, the Research Corporation, and the Texas Advanced Research Program.


format: text

text:
# Simulation of the AGILE Gamma-Ray Imaging Detector Performance: Part I

## 1 Introduction

AGILE is an Italian Space Agency (ASI) Small Scientific Mission dedicated to high-energy astrophysics . The AGILE instrument is designed to detect and image photons in the $`30`$ MeV–50 GeV and 10–40 keV energy bands, with excellent spatial resolution and timing capability, and a very large field of view covering $``$ 1/5 of the entire sky at energies above 30 MeV . Primary scientific goals include the study of AGNs, gamma-ray bursts, Galactic sources, unidentified gamma-ray sources, diffuse Galactic and extragalactic gamma-ray emission, high-precision timing studies and Quantum Gravity testing .

The AGILE Gamma-Ray Imaging Detector (GRID) consists of a Silicon-Tungsten Tracker, a Cesium Iodide Mini-Calorimeter (MC), an Anticoincidence system (AC) made of segmented plastic scintillators, and a Data Handling system (DH). The GRID is sensitive in the energy range $`30`$ MeV–50 GeV, and is designed to achieve an optimal angular resolution (source location accuracy $``$ 5’–20’ for intense sources), a very large field-of-view ($`3`$ sr), and a sensitivity comparable to that of EGRET for on-axis (and substantially better for off-axis) point sources. The optimal performance of the AGILE GRID requires an efficient on-board selection of events aimed at maximizing the background rejection and the transmission of cosmic photon events. This task is carried out by the Data Handling (DH) System.

In the current (Paper I) and companion  papers we present the main results concerning the simulation of GRID events induced by the particle/albedo-photon background and by the cosmic gamma-rays. We study in detail different levels of the GRID on-board trigger and data processing. For the Montecarlo simulations we used the numerical code GEANT (version 3.21)  and related subroutines. Our numerical results are in very good agreement with CERN beam-test results obtained by our group with electrons and photons in the 0.1–1 GeV energy range .

## 2 Cosmic gamma-ray signal and expected particle and albedo-photon background

### 2.1 Cosmic gamma-ray flux

The cosmic photon flux on the AGILE-GRID detector will be dominated by the diffuse gamma-ray emission, typically larger near the Galactic plane by a factor of 10 with respect to the high-Galactic-latitude. This latter emission has a flux of $`\mathrm{\Phi }10^5\mathrm{photons}/\mathrm{cm}^2/\mathrm{s}/\mathrm{sr}`$ for photons above 100 MeV. Expected gamma-ray event rates for the AGILE-GRID depend on the pointed fields and can be given in the ranges (dominated by the diffuse gamma-ray emission): 0.02 - 0.1 photons/s for energies larger than 100 MeV and 0.1 - 1 photons/s for energies between 10 and 100 MeV. The contribution of point-like gamma-ray sources adds to these rates without changing their order of magnitude.

In the simulations we considered 4 classes of cosmic gamma-ray photons, divided in low-energy photons (25–50 MeV) and high-energy photons (400–1000 MeV); we used a power-law energy spectrum of index -2, and run simulations for two extreme cases: an input geometry of incidence angle between 0 and 10 degrees off-axis, and an input geometry of incidence angle between 50 and 60 degrees off-axis.

### 2.2 Charged particle background

An equatorial orbit near 550 km of altitude for the AGILE mission will provide a relatively low-background environment. The charged particle background for this orbit is known to be relatively stable, with an increase by a factor 10-100 near the South Atlantic Anomaly (typically influencing about 10% of the orbital duration). Taking into account data from SAS-2  and Beppo-SAX  missions, out of the South Atlantic Anomaly we expect an average rate of charged particle background above $`1`$ MeV of $`0.3`$ particles cm-2 s-1 (mostly electrons/positrons with a $``$10% contribution of protons).

In our simulations, based on spectral data from AMS  and MARYA , we considered the following background components: electrons (ELE), positrons (POS), primary high-energy protons (HE PROT) and secondary low-energy protons (LE PROT). We made the following assumptions: 1) we considered an “average equatorial position” and we assumed the AGILE z-axis pointed towards the sky for an unocculted GRID field of view; 2) we considered a spectrum averaged over the magnetic latitudes $`|\lambda _{mag}|=011^{}`$ for secondary protons, and $`|\lambda _{mag}|=017^{}`$ for electrons and positrons (values compatible with the AGILE equatorial orbit), regardless of longitude dependence; 3) we considered an isotropic distribution on a virtual sphere containing the satellite for particles below the geomagnetic cutoff (albedo leptons and secondary protons), and an isotropic distribution emitted from the upper hemisphere for the primary protons above the cutoff. The adopted particle spectra, which are fits of the AMS and MARYA data (calculated in ) are shown in Fig. 1. Their fluxes are normalized to the values indicated in Tab. 1 and the total incidence rate on a surface is assumed to be $``$ 0.3 particles $`\mathrm{cm}^2\mathrm{s}^1`$ for particle kinetic energies above $`10`$ MeV.

### 2.3 Albedo gamma-ray background

The interaction of the charged cosmic-rays with the upper atmosphere induces a relatively strong and non-isotropic gamma-ray background. It peaks near the Earth horizon (corresponding to a zenith angle $`\theta =112^{}`$ for an orbit of 550 km altitude), and has a characteristic East-West asymmetry (by a factor of $`4`$ in intensity near 40 MeV) .

In our simulations, we used a simplified but representative model of the albedo-photon emission, and we considered two geometries: (1) albedo photons reaching the GRID ”from below”, a situation valid for the unocculted portions of the AGILE orbit (ALB-1 PHOT); (2) Earth covering half of the AGILE FOV, with a large fraction of albedo photons reaching the GRID ”from one lateral side” (ALB-2 PHOT). Based on balloon flight data  and SAS-2 data , we have assumed the following simplified spectrum representing the average of the emission over the whole solid angle subtended by the Earth surface ($`4`$ steradians) at the height of 550 km:

$`{\displaystyle \frac{dN}{dSdtdEd\mathrm{\Omega }}}=\{\begin{array}{cc}0.08E^{1.4}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{MeV}^1\hfill & \text{for }1\mathrm{MeV}<\mathrm{E}<10\mathrm{MeV}\hfill \\ 0.5E^{2.2}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{MeV}^1\hfill & \text{for }10\mathrm{MeV}<\mathrm{E}<100\mathrm{GeV}\hfill \end{array}`$ (3)
with a total flux $`=0.15\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1`$ in the energy range 1 MeV - 100 GeV.

### 2.4 Background rejection

The ratio of charged particle events (penetrating the AC) to cosmic gamma-ray photon events is typically of order $`10^310^4`$ for photons of energy 20-100 MeV, and $`10^410^5`$ for photons above 100 MeV. Albedo photons are larger than cosmic gamma-ray photons by factors between 10 and 100 (depending on the pointing geometry and energy range).

Clearly, is necessary to ‘filter’ the GRID events. This task is carried out by a hardware-implemented fast logic (“Level-1 trigger”), and by a set of asynchronous software algorithms and CPU processing (”Level-2 processing”). The adopted requirements for the DH processing depend on the downlink telemetry rate, assumed to be 512 kbit/s during the 10 minute duration of satellite visibility from the ASI ground station in Malindi (Kenya). A main aim of the AGILE DH system is to provide an on-board filtering of events reducing the background rate to an acceptable value within a factor of 10–100 of the cosmic gamma-ray photon rate. A final filtering of the particle and albedo-photon background will be carried out on the ground by a dedicated software. The charged particle and the albedo-photon background models described in the previous sections were used in Montecarlo simulations in order to optimize the DH processing.

## 3 The AGILE Model and GRID detection assumptions

### 3.1 The AGILE Model

A correct mathematical and physical model of the AGILE payload and spacecraft with all their relevant components is crucial to evaluate the GRID behaviour and performance. A general description of the AGILE instrument can be found in . Simulation programs consider the following AGILE components (see Fig. 2): the spacecraft (bus) MITA and the AGILE payload, consisting of the CsI Mini-Calorimeter, the Silicon-Tungsten Tracker, the Anticoincidence system, the X-ray detector (Super-AGILE), the thermal shield, the Mechanical structure and the lateral electronics boards. The AGILE components were modeled as follows:

* MITA spacecraft
The spacecraft (bus) MITA is represented by a box, made of Carbon-fibers, containing 10 equidistant Aluminum layers inserted in thin (100 $`\mu \mathrm{m}`$) Carbon-fiber structures; every Aluminum layer is 1.1 cm thick for a total weight of 128 kg.
* Mini-Calorimeter
The Mini-Calorimeter is modeled by two planes, each containing 16 CsI bars, oriented as the y-axis in the upper plane and as the x-axis in the lower plane. The CsI bars are 1.4 cm thick and have a pitch of 2.4 cm. They are inserted in a Carbon-fiber structure 0.1 cm apart one each other. The distance between the X layer and the Y layer is assumed to be 0.1 cm. The Mini-Calorimeter is placed directly above the bus MITA, at a 0.5 cm distance from the Silicon Tracker.
* Silicon Tracker
The Silicon Tracker is made of 14 detection planes, with a distance between consecutive planes of 1.6 cm and a distance between the X and Y layer in each plane of 0.16 cm; the first 12 planes contain also a Tungsten layer each ($`245\mu \mathrm{m}`$ thick). Each Silicon layer is composed of 16 tiles of $`9.5\times 9.5\mathrm{cm}^2`$ each, with a pitch of $`121\mu \mathrm{m}`$ (each tile contains 768 strips) and a thickness of $`410\mu \mathrm{m}`$. In our model there are also the honeycomb support, the Aluminum inserts connecting the Tracker planes with the main structure of the satellite, and the Front-End electronics chips (3 TAA1 for each ladder). Every Tracker plane is composed by a Carbon-fiber structure ($`500\mu \mathrm{m}`$ thick) supporting an Aluminum honeycomb. Under the Carbon-fiber are glued the tungsten layer (for the first 12 planes), one layer of Kapton ($`100\mu \mathrm{m}`$ thick) with a Copper coating ($`15\mu \mathrm{m}`$ thick) and the X Silicon layer. Under the X layer there are the Y layer and another Kapton-Copper layer, which in the reality are glued to the next plane structure. On the same structure are glued also the ceramic hybrids which sustain the TAA1 Silicon chips. The global structure is composed by 14 X and Y detectors layers and by 15 structure supports; the total Tracker height is 24 cm.
* Anticoincidence system
The AC system is made of a top panel of plastic scintillator ($`54\times 54\times 0.5\mathrm{cm}^3`$) and 3 overlapping panels for each side of the AGILE Tracker (0.6 cm thick, 18.1 cm large and 44.4 cm high). The AC plastic panels are inserted in a Carbon-fiber structure ($`500\mu \mathrm{m}`$ thick) and are supported by some Aluminum inserts; the trapezoidal form of the lateral panels is described in a proper way. We have inserted also a schematic description of the photomultipliers and their casings. The AC system surrounds the Mini-Calorimeter, the Silicon Tracker and Super-AGILE, and is positioned above the bus MITA. Fig. 3 shows the adopted model.
* X-ray detector (Super-AGILE)
The X-ray detector is made substantially by these elements: a Silicon detection layer, a collimation system, a gold mask. A Tungsten ring is positioned below the Si detection plane in order to reduce the diffuse photon background. The detection layer is divided in 16 tiles with the same characteristics of the Tracker tiles, described in a realistic way, with the Silicon microstrips and the FrontEnd electronics. We have modeled also the honeycomb and Carbon-fiber support structure; the detectors and the support of the coded mask are glued on this structure. The Silicon detector is placed over the Tracker, while the distance between the Silicon detector and the gold mask ($`90\mu \mathrm{m}`$ thick) is fixed to be 14 cm. Also the collimation system, made of an ultra-light Carbon-fiber structure ($`500\mu \mathrm{m}`$ thick), is described in great detail (Fig. 4). Its internal panels are coated with a $`75\mu \mathrm{m}`$ thick gold layer, in order to reduce the contribution from the diffuse cosmic X-ray background. The distance between the gold mask and the top AC panel is 0.5 cm.
* Thermal shield
The thermal shield is modeled as an external Teflon layer (1 mm thick) plus an internal Aluminum layer (10 $`\mu \mathrm{m}`$ thick). It covers the Mini-Calorimeter, the Tracker, Super-AGILE and the AC system.
* Mechanical Structure
The description of the payload is completed by a preliminary model for the general support structure: four Aluminum legs connected by other Aluminum bars.
There is also a preliminary description of the Front End Electronics (FEE) and of the trigger vertical boards, one for each lateral side of the payload, and one for the DH box, under the Mini-Calorimeter. They are described as some layers of Silicon, ceramics and Copper.

### 3.2 Silicon Tracker capacitive coupling, floating strip readout, and Gaussian noise

The electronic system for the AGILE Silicon Tracker is based on the “floating strip” readout. This implies that only one out of two contiguous Si-microstrips (each of 121 $`\mu \mathrm{m}`$ size) is read by dedicated electronic devices (TAA1 chips). However, the capacitive coupling between contiguous strips allows not to loose any information. In the Montecarlo code the capacitive coupling is simulated calculating the total energy release in every readout strip by the weighted sum as obtained in Ref. . For the $`k`$-th readout strip, the total energy release taking into account neighbour strips is $`E_s(k)=E_i+0.38(E_{i1}+E_{i+1})+0.115(E_{i2}+E_{i+2})+0.095(E_{i3}+E_{i+3})+0.045(E_{i4}+E_{i+4})+0.035(E_{i5}+E_{i+5})`$, where $`E_i`$ is the energy release in the $`i`$-th strip that can be either readout or floating, and $`k=1+(i1)2`$. This scheme for the capacitive coupling simulation is in agreement with test-beam results obtained from our group for a variety of incidence angles .

Instrumental noise was simulated as a Gaussian distribution with $`\sigma =5`$ keV. We found that for a detection threshold of 1/4 MIP$`27`$keV it has no macroscopic effects on the trigger rates. All the simulation results presented in this document have been obtained using the geometry model described above, taking into account the capacitive coupling and the floating strip readout, but without considering the noise that has negligible effects on the results.

## 4 Level-1 trigger

The Level-1 trigger for the Agile GRID is a hardware-implemented fast logic required to be very rapid ($``$ 2 $`\mu \mathrm{s}`$), conceptually simple, easily implementable by a dedicated hardware, and using parameters reconfigurable by Telecommands. With Montecarlo simulations we tested different trigger configurations, using the event classes (cosmic gamma-ray photons, background particles and Earth albedo photons) described in section 2. Here we present the main results.

Our best strategy for the Level-1 logic is based on the combined use of signals from AC panels and of the ratio R, defined as the ratio between the total number of fired TAA1 and the total number of fired X and Y views: $`R=(\text{total n. of fired TAA1})/(\text{total n. of fired X and Y views})`$. In what follows, we call ”$`R`$–trigger” the logic (implemented by a dedicated electronic component) regulating the use of both AC panels and $`R`$ threshold values.

### 4.1 Level-1 trigger configurations

Here is a short description of the different Level-1 trigger configurations that we studied:

* PLA = events which give hits in at least 3 out of 4 consecutive planes (X OR Y view)
* TOP = PLA events which pass the TOP AC veto
* LAT = TOP events with signals in 0 or 1 lateral AC panels, and TOP events which give signals in 2 consecutive AC panels or in 2 AC panels on the same side
* LSI = LAT events with signals in 0 lateral AC panels, and LAT events with signals in 1 or 2 AC panels but with NO signal in the last silicon plane
* R11G = LAT events with signals in 0 lateral AC panels, and LAT events with signals in 1 or 2 AC panels and $`R>1.1`$

PLA, TOP and LAT are consecutive steps, while LSI and R11G are the two final alternative trigger configurations. For comparison, in Tables 2, 3 and 4 we reported also the number of events characterized by primary particles or photons reaching the Tracker volume (TRA), a box of $`38.06\times 38.06\times 21.078`$ $`\mathrm{cm}^3`$ which includes the Tracker planes from the top sheet of the first tungsten layer to the bottom sheet of the last Silicon-y plane.

In the case of gamma-ray photons it is important to observe that:

* gamma-photons giving a signal in the Tracker may convert in the Tracker volume, or in the Super-AGILE volume (mask or collimators). They can also convert elsewhere (mechanical structure, Mini-Calorimeter or bus MITA), or pass through the instrument without converting. We define as “good photons”, those for which there is high probability to reconstruct the incident direction. By definition they are “Tracker-converted photons”;
* because of the presence of the Super-AGILE structure, some gamma-photons potentially able to enter into the Si-Tracker fiducial volume, because of their incidence angles and energies, can convert in the Super-AGILE structure. Their secondaries (electron and positron) may never arrive into the Tracker volume. On the other hand, it can happen that some photons, that geometrically would not enter into the Tracker volume, convert in the Super-AGILE structure, and their secondaries give a signal in the Tracker then inducing a GRID background.

It is then useful to indicate with the suffix “TC” the photons converted in the Tracker volume, and with the suffix “SC” the photons converted in the Super-AGILE structure. In the case of gamma-ray photons we also distinguish: a) photons that theoretically, from a geometrical point of view, might enter into the Tracker volume (TRA-TH), b) primary photons actually entered in the Tracker volume (PTRA), c) primary photons converted in the Tracker volume (PTCON).

### 4.2 Results and discussion

Tabs. 2, 3 and 4 and Fig. 5 summarize our results in terms of alternative or progressive event selections. The TOP AC veto is useful to reject downgoing charged particles; the LAT AC veto discriminates events giving a signal either entering than exiting the detector (charged particles) from events which give signals only when they escape from the Tracker (photons into charged pairs). The R-trigger is a method optimized to reject non-interacting charged particles, based on a simplified idealization: when a non-interacting particle enters in the Tracker, only one track is expected; if the detector is crossed by a photon creating an electron-positron pair, two tracks are expected. Simulation results show that the R-trigger works very well with protons, and gives acceptable results also for electrons and positrons. In general, the method succeeds in discriminating photons from charged background particles very efficiently, sometimes better than the AC system. Clearly, the use of the $`R`$-trigger does not reject protons generating many secondaries or e+/e- producing secondaries and electromagnetic showers. A trigger configuration based only on the AC system does not reject low energy protons stopping inside the Tracker. The two methods are based on different “philosophies” and they both give good results in different situations.

From our simulations it appears that the best choice is a combination of the two strategies and that the most efficient Level-1 trigger configuration is R11G, rejecting $`93\%`$ of background particles without affecting significantly the gamma-ray detection rate.

## 5 Level-1.5 Data Processing

The additional “Level-1.5” data processing uses the Tracker discretization in terms of TAA1 chips ($`3.1\times 3.1`$ $`\mathrm{cm}^2`$ square elements), together with the complete AC and MC information. We investigated the efficiency of alternative Level-1.5 Trigger options in rejecting background charged particles and here we present the main results.

### 5.1 Level 1.5 trigger configurations

We considered the following options:

* DIS =The DIS option is a simplified track reconstruction and is based on computing the distance $`D`$ of the fired TAA1s from the fired AC lateral panel. The basic idea is that a charged particle is expected to have increasing $`D`$ values as a function of increasing plane, while a gamma-ray photon is expected to have a decreasing $`D`$ function (in the restricting hypothesis that they both come downward and in the absence of strong “hard” scattering of photon-created pairs). The parameter DIS is defined as: $`DIS=D_{fp}D_{lp}`$ where $`D_{fp}`$ is the distance of the closest fired TAA1 to the fired AC lateral panel in the first plane, while $`D_{lp}`$ is the distance of the closest fired TAA1 to the fired AC lateral panel in the last plane. We require $`DIS0`$ for good events. This option is applied only if there are fired AC lateral panels.
* RUD10 (or RUD12) = The fired Tracker views are divided in two groups (UP and DOWN) and the ratio R is computed separately for the two groups ($`R^{up}`$ and $`R^{down}`$). The ratio $`RUD=R^{up}/R^{down}`$ can help in discriminating good gamma-photons, that come downward, from secondary photons generated by charged particles, that come from the bottom toward the top of the detector. With the option RUD10 are rejected all the events with $`RUD>1.0`$, while with the option RUD12 are rejected all the events with $`RUD>1.2`$.
* VFORM = This option, as RUD10 and RUD12, can help in discriminating primary photons going downward from secondary photons coming upward. It computes the separation $`S`$ between fired TAA1 on the same view, divides the fired planes in two groups (UP and DOWN) and compares the mean separation of the UP planes with the mean separation of the DOWN planes. The parameter $`V`$ is defined as $`V=S_{down}S_{up}`$ and events with $`V<0`$ are rejected.
* CGAP = This option rejects events with a gap of more than 2 planes (4 views) between the MC and the last fired Tracker view. It is applied only if there is a signal from the MC and it is supposed to reject events due to $`e^+`$ or $`e^{}`$ which stop in the MC and generate a secondary photon which enters the Tracker upward and creates a couple only after passing 4 or more views.
* TGAP = This option rejects events with a gap of more than 2 planes (4 views) in the Tracker. The parameter T\_GAP is defined, for a single event, as the maximum number of no-fired views between two non consecutive fired views.

### 5.2 Results and Discussion on the Level-1.5 Trigger Processing

Tables 6 and 6 and Fig. 6 show the simulation results. The simple DIS algorithm is the most efficient in rejecting background particles without losing too many good photons. This algorithm is applied only if there is at least one fired AC lateral panel; most of background charged particles verify this condition, while good photons often do not hit any lateral AC panel.

Options RUD10, RUD12 and VFORM are not very efficient because the distributions of the RUD and V parameters are not significantly different for background particles and for good photons. The C\_GAP option is not efficient because the majority of background electrons, positrons and albedo protons do not give a signal in the Mini-Calorimeter. There are no significant differences between the background particle T\_GAP distributions and the photon T\_GAP distribution, probably because this parameter is influenced by secondary particles.

## 6 Discussion and Conclusions

By using a complete and detailed model of the AGILE instrument, we simulated with GEANT the performance of the on-board Level-1 trigger logic. We obtain satisfactory results from our optimization study of the Level-1 and 1.5 trigger logic for the cosmic gamma-ray signal (with photon energies between $`30`$ MeV and 50 GeV) and for the rejection of the particle and albedo-photon backgrounds in a 550 km equatorial orbit. The optimal Level-1 (and therefore necessarily simple) event selection algorithms are the R11G and DIS procedures that can be easily hardware-implemented and are very fast ($`20\mu `$s). Table 7 summarizes the main results of our paper: the total background rates expected for two different Earth occultation geometries. At the Level-1$`+`$1.5 the background component completely dominates the event rate over the cosmic gamma-rays by a factor between 10 and 100, and the total background event rate is near 100 $`\mathrm{s}^1`$, a value which can be sustained by the on-board AGILE Data Handling . Therefore the Level-1 logic presented here is adequate for the AGILE GRID.

We also note that the albedo photons contribute a significant fraction (between $`30`$% and more than $`50`$% depending on the Earth’s position with respect with the GRID) of the total background rate. Since these photons are dominated by 10–100 MeV gamma-rays, Level-1 techniques are not adequate in substantially reducing their event rate, and a special Level-2 data processing is required,which will be the subject of a forthcoming paper.

Our results are of general interest for space detectors similar to the AGILE GRID. Two conclusions are of particular interest: (1) the fact that low-energy leptons (magnetospheric trapped electrons and positrons) dominate the particle background requires special attention and a proper event selection logic; (2) Earth albedo-photons constitute a substantial fraction of the residual events passing the Level-1 trigger stage.

## 7 Acknowledgments

Results presented in this paper are based on joint work with members of the AGILE Team. In particular, we warmly thank G. Barbiellini, P. Picozza and A. Morselli for special support, M. Prest, E. Vallazza and G. Fedel for exchange of information on the Silicon Tracker performance and experimental data, and the AGILE Simulation and Theory Group for many discussions. We also thank A. Pellizzoni and P. Lipari for collaboration on several background source function algorithms used in our simulations.

The current work was carried out at the University of Rome ”Tor Vergata”, University of Ferrara and CNR and INFN laboratories, under the auspices and partial support of the Agenzia Spaziale Italiana.

FIGURE CAPTIONS:
Fig.1: Charged particle background for an equatorial orbit of 550 km (results from ref. ) that we assumed in our calculations. Electron, positron and proton spectra were obtained by combined AMS-Shuttle Flight data and MARYA data. All particle components have an approximately isotropic distribution of incidence angles on the GRID, except the primary high-energy protons above the geomagnetic cutoff of $`6`$ GeV.
Fig.2: AGILE payload: from the bottom the different volumes are: the DH box, the 2 planes of the Mini-Calorimeter, the 14 layers of the Silicon Tracker, Super-AGILE, the top AC panel. Also the lateral AC panels with the photomultipliers and the thermal shield are visible. The overall dimensions are: length of 62.55 cm, width of 62.55 cm and height of 54.02 cm. The GRID active volume dimensions are: $`38.06\times 38.06\times 21.078`$ $`\mathrm{cm}^3`$ including the Tracker planes from the top sheet of the first Tungsten layer to the bottom sheet of the last Silicon-y plane.
Fig.3: Anticoincidence system. The overall dimensions are: lenght of 62.25 cm; width of 62.25 cm; height of 53.87 cm.
Fig.4: Super-AGILE. Overall dimensions: lenght of 44 cm; width of 44 cm; height of 14.77 cm.
Fig.5: Particle background rates for different Level-1 trigger configurations and Level-1 trigger efficiency for different photon classes. Event selection based on TRA, PLA, TOP and LAT are sequentially applied to all events. Further selections based on LSI or R11G are instead mutually exclusive. The suffix “TC” denotes “Tracker-converted”.
Fig.6: These plots show the efficiency of the Level-1.5 Trigger alternative options in rejecting particle background and in detecting gamma-photons (DISRUD means (DIS+RUD\_10) and DISCG means (DIS+C\_GAP)).


text:
# Two-dimensional negative donors in magnetic fields

## I Introduction

The properties of $`D^{}`$ centers in narrow quantum wells have become a subject of considerable interest during the past years Phelps ; Larsen1 ; Larsen2 ; XuWang ; Sandler ; TaoPang ; Louie ; Zhu1 ; Zhu2 . In very narrow wells this system can be considered as a two-dimensional counterpart of the atomic $`\mathrm{H}^{}`$ ion. A frequently used compound for current experimental investigations (see for example ref.Jiang ) of such systems are layers of GaAs/AlGaAs. The high mobility i.e. small effective mass of the electrons and the comparatively large dielectricity constant of this semiconductor material allow us to study strong magnetic field effects in the laboratory. Both the two-dimensional character of the motion of the electrons and the external magnetic field makes the ground state of the system more tightly bound than that of the field-free threedimensional $`\mathrm{H}^{}`$ ion. Considering neutral donors $`D^0`$ the planar electron density of the 2d donor is enhanced significantly compared to that of the 3d donor. The form of the corresponding charge distribution makes the interaction of an additional distant charge with the neutral donor very different for the three-dimensional compared to the two-dimensional situation. In three dimensions it is well-known that, for a fixed center (nucleus), the combination of a magnetic field with the 3d long-range attractive polarization forces lead to an infinite number of bound states for the negative donor Avr81 ; Bez00 ; Bez01 ; AlH00 . In two dimensions the long-range interaction is of repulsive character. Indeed, it was shown by Larsen and McCann Larsen1 , that an electron situated far from a 2d $`D^0`$ center experiences an overall repulsive potential which prevails both the attractive Coulomb attraction due to the center and the attractive polarization forces. For distances of the order of the extension of the neutral donor, exchange and correlation effects play a major role with respect to the binding mechanism of the additional electron.

Two-dimensional $`D^{}`$ centers in magnetic fields are described by a three-dimensional Schrödinger equation. The finite difference numerical method Ivanov86 ; Ivanov85 ; Ivanov88 ; Ivanov94 ; Ivanov98 ; IvaSchm2001b ; IvaSchm2001c ; Ivanov86a  allows us to solve this Schrödinger equation without any simplifications or approximations with respect to the geometry of the wave function, correlation effects etc. Due to this property of our computational approach it will be here possible to obtain a definite answer on the ground state energy and the number of bound states of the two-dimensional $`D^{}`$ center in magnetic fields. The critical values of the field strengths that are associated with the appearance of the bound states will be determined. One of the major motivations for the present work are the different theoretical values for the ground state energy of a two-dimensional negative donor with and without magnetic field existing in the literature. Variational calculations Phelps ; Larsen2  on the onehand and Monte Carlo simulations Louie  on the otherhand yield uncompatible predictions for the ground state energies. This puzzle will be resolved within the present work.

## II The Schrödinger equation and the method of solution

The Hamiltonian of our two-dimensional system of two interacting electrons with an effective mass $`m`$ and a singly charged positive ion placed into a magnetic field perpendicular to the plane can be written in Cartesian coordinates $`(x,y)`$ as

$`H={\displaystyle \frac{1}{2}}(_1^2+_2^2)i{\displaystyle \frac{\gamma }{2}}\left(y_1{\displaystyle \frac{}{x_1}}+x_1{\displaystyle \frac{}{y_1}}y_2{\displaystyle \frac{}{x_2}}+x_2{\displaystyle \frac{}{y_2}}\right)+{\displaystyle \frac{\gamma ^2}{8}}(r_1^2+r_2^2)`$
$`{\displaystyle \frac{1}{r_1}}{\displaystyle \frac{1}{r_2}}+{\displaystyle \frac{1}{r_{12}}}`$ (1)
Here we use the units $`a_{\mathrm{eff}}=\mathrm{}^2ϵe^2m^1=ϵm^15.29\times 10^9`$cm for the length, $`E_{\mathrm{eff}}=e^4\mathrm{}^2mϵ^2=m^{}ϵ^227.2`$eV for the energy and $`B_{\mathrm{eff}}=ce^3m^2\mathrm{}^3ϵ^2=m^2ϵ^22.3505\times 10^5`$T for the magnetic field strength. $`\gamma =B/B_{\mathrm{eff}}`$, $`m_0`$ is the mass of the free electron, $`m^{}=m/m_0`$, $`ϵ`$ is the dielectricity constant of the semiconductor material, $`r_1=(x_1^2+y_1^2)^{1/2}`$, $`r_2=(x_2^2+y_2^2)^{1/2}`$, and $`r_{12}=\left|\stackrel{}{r}_1\stackrel{}{r}_2\right|`$. The spin terms are omitted.

In contrast to variational calculations we are not biased by using a particular ansatz for the variational wave function but will perform a full grid solution to the Schrödinger equation that allows to control and estimate the remaining minor deviation from the exact eigenfunctions and eigenvalues. In the following we provide an outline of our approach which embodies various properties and transformations of the Hamiltonian and the corresponding Schrödinger equation.

Transforming the electronic degrees of freedom to polar coordinates $`(r_1,\varphi _1)`$ and $`(r_2,\varphi _2)`$ yields for the Hamiltonian (1)

$`H={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{^2}{r_1^2}}+{\displaystyle \frac{1}{r_1}}{\displaystyle \frac{}{r_1}}+{\displaystyle \frac{^2}{r_2^2}}+{\displaystyle \frac{1}{r_2}}{\displaystyle \frac{}{r_2}}\right){\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{r_1^2}}{\displaystyle \frac{^2}{\varphi _1^2}}+{\displaystyle \frac{1}{r_2^2}}{\displaystyle \frac{^2}{\varphi _2^2}}\right)i{\displaystyle \frac{\gamma }{2}}\left({\displaystyle \frac{}{\varphi _1}}+{\displaystyle \frac{}{\varphi _2}}\right)`$
$`+{\displaystyle \frac{\gamma ^2}{8}}(r_1^2+r_2^2){\displaystyle \frac{1}{r_1}}{\displaystyle \frac{1}{r_2}}+{\displaystyle \frac{1}{(r_1^2+r_2^22r_1r_2\mathrm{cos}(\varphi _1\varphi _2))^{1/2}}}`$ (2)
The eigenfunctions of this Hamiltonian are the eigenfunctions of the Hamiltonian (1) if we require them to be periodic with period $`2\pi `$ with respect to the variables $`\varphi _1`$ and $`\varphi _2`$. The eigenfunctions of the Hamiltonian (2) should also be eigenfunctions of the $`z`$-projection of the total orbital angular momentum operator which is a conserved quantity ie. $`𝐥_z\mathrm{\Psi }=M\mathrm{\Psi }`$ with

$`𝐥_z=i{\displaystyle \frac{}{\varphi _1}}i{\displaystyle \frac{}{\varphi _2}}`$ (3)
Introducing the coordinates $`r_1,r_2,\varphi ,\mathrm{\Phi }`$ where

$`\varphi =\varphi _1\varphi _2`$ (4)
and

$`\mathrm{\Phi }={\displaystyle \frac{\varphi _1+\varphi _2}{2}}`$ (5)
we have

$`𝐥_z=i{\displaystyle \frac{}{\mathrm{\Phi }}}`$ (6)
For these coordinates the Hamiltonian takes on the appearance

$`H={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{^2}{r_1^2}}+{\displaystyle \frac{1}{r_1}}{\displaystyle \frac{}{r_1}}+{\displaystyle \frac{^2}{r_2^2}}+{\displaystyle \frac{1}{r_2}}{\displaystyle \frac{}{r_2}}\right){\displaystyle \frac{1}{8}}\left({\displaystyle \frac{1}{r_1^2}}+{\displaystyle \frac{1}{r_2^2}}\right){\displaystyle \frac{^2}{\mathrm{\Phi }^2}}`$
$`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{r_1^2}}{\displaystyle \frac{1}{r_2^2}}\right){\displaystyle \frac{^2}{\mathrm{\Phi }\varphi }}{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{r_1^2}}+{\displaystyle \frac{1}{r_2^2}}\right){\displaystyle \frac{^2}{\varphi ^2}}i{\displaystyle \frac{\gamma }{2}}{\displaystyle \frac{}{\mathrm{\Phi }}}`$
$`+{\displaystyle \frac{\gamma ^2}{8}}(r_1^2+r_2^2){\displaystyle \frac{1}{r_1}}{\displaystyle \frac{1}{r_2}}+{\displaystyle \frac{1}{(r_1^2+r_2^22r_1r_2\mathrm{cos}\varphi )^{1/2}}}`$ (7)
Eigenfunctions of the operator $`𝐥_z`$ read as follows

$`\mathrm{\Psi }(r_1,r_2,\varphi ,\mathrm{\Phi })=e^{iM\mathrm{\Phi }}\psi (r_1,r_2,\varphi )`$ (8)
We therefore have

$`H\mathrm{\Psi }=e^{iM\mathrm{\Phi }}h\psi `$ (9)
with

$`h={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{^2}{r_1^2}}+{\displaystyle \frac{1}{r_1}}{\displaystyle \frac{}{r_1}}+{\displaystyle \frac{^2}{r_2^2}}+{\displaystyle \frac{1}{r_2}}{\displaystyle \frac{}{r_2}}\right)+{\displaystyle \frac{M^2}{8}}\left({\displaystyle \frac{1}{r_1^2}}+{\displaystyle \frac{1}{r_2^2}}\right)`$
$`{\displaystyle \frac{iM}{2}}\left({\displaystyle \frac{1}{r_1^2}}{\displaystyle \frac{1}{r_2^2}}\right){\displaystyle \frac{}{\varphi }}{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{r_1^2}}+{\displaystyle \frac{1}{r_2^2}}\right){\displaystyle \frac{^2}{\varphi ^2}}+{\displaystyle \frac{M\gamma }{2}}`$
$`+{\displaystyle \frac{\gamma ^2}{8}}(r_1^2+r_2^2){\displaystyle \frac{1}{r_1}}{\displaystyle \frac{1}{r_2}}+{\displaystyle \frac{1}{(r_1^2+r_2^22r_1r_2\mathrm{cos}\varphi )^{1/2}}}`$ (10)
This Hamiltonian depends on three degrees of freedom $`(r_1,r_2,\varphi )`$ and the corresponding Schrödinger equation can therefore be solved for the eigenfunctions and eigenvalues by applying our finite-difference method Ivanov86 ; Ivanov85 ; Ivanov88 . The coordinates $`r_1`$ and $`r_2`$ are in the domain

$`0<r_1<+\mathrm{},0<r_2<+\mathrm{}`$ (11)
with zero boundary conditions for infinity. The corresponding domain and boundary conditions for the angle $`\varphi `$ have to be analyzed in more detail.

From equation (5) one can conclude that solutions of the Schrödinger equation belonging to the Hamiltonian (10) should be considered in the domain $`2\pi \varphi 2\pi `$ with cyclic boundary conditions. This leads to eigenfunctions that can be either symmetric or antisymmetric with respect to interparticle exchange. The first correspond to spin singlet states whereas the second yield spin triplet states. Beyond this we encounter additional eigenfunctions of the operator (10) that are not eigenfunctions of our initial Hamiltonian (1) i.e. these eigenfunctions do not describe physical solutions. Indeed, exploiting certain symmetries of the wavefunction it is possible to confine the domain of the angle $`\varphi `$ to $`0\varphi 2\pi `$ for the Schrödinger equation belonging to the Hamiltonian (10). As a consequence the non-physical solutions are excluded. For the states with $`M=0`$ solutions can be obtained in an even smaller domain $`0\varphi \pi `$.

The numerical finite-difference method employed here to solve the eigenvalue problem for the Hamiltonian (10) is a modification of the approach developed in previous works devoted to atoms in strong magnetic Ivanov88 ; Ivanov94 ; IvaSchm2001b  and electric fields (seeIvanov98 ; IvaSchm2001c  and references therein) and contains, in particular, technical aspects used in three-dimensional problems Ivanov85 ; Ivanov86a . Our computational procedure consists of the following main steps. The nodes have to be chosen in the spatial domain $`\mathrm{\Omega }:0r_1<+\mathrm{}`$, $`0r_2<+\mathrm{}`$, $`0\varphi 2\pi `$ (or $`0\varphi \pi `$). The values of the wavefunctions at the positions of the nodes are the numerical representation of the solutions of the initial differential equations. The nodes of the three-dimensional mesh in the space $`(r_1,r_2,\varphi )`$ are placed at all the points with coordinates $`(r_{1i},r_{2j},\varphi _k)`$ where $`r_{1i}`$, $`r_{2j}`$, and $`\varphi _k`$ are sets of the mesh node coordinates along the corresponding axes. For the coordinate $`\varphi `$ it is natural to use a uniform mesh with nodes at $`\varphi _k=2\pi (k1/2)/N`$ for the region $`0\varphi 2\pi `$ (case $`M0`$) and $`\varphi _k=\pi (k1/2)/N`$ for the region $`0\varphi \pi `$ (case $`M=0`$), where $`N`$ is the number of nodes in the direction $`\varphi `$. For $`r_1`$ and $`r_2`$ we have used non-uniform distributions of nodes similar to those described in ref.IvaSchm2001b , which cover the infinite domains $`0r_1<+\mathrm{}`$, $`0r_2<+\mathrm{}`$ as $`N\mathrm{}`$.

We employ the inverse iteration method to obtain the eigenfunctions and energy eigenvalues. This requires solving a system of linear equations with a matrix that is a finite-difference approximation to the Hamiltonian. The solution of these equations is particularly simple if the matrix has a block-tridiagonal structure. The latter can be achieved using the simplest three-point approximation for the derivatives for one of the coordinates. The choice of this coordinate is dictated by obtaining a convenient form of representation for the boundary conditions. For the two other coordinates we are free to apply more precise higher order approximations to their derivatives. The final values for the energy (and other observables) are provided by using the Richardson extrapolation technique for the corresponding results emerging from a series of geometrically similar meshes with different number of nodes. Using this approach we achieve a major increase of the numerical precision and, in particular, we obtain together with each numerical value a reliable estimate of its precisionIvanov86 ; Ivanov88 . Typically, meshes used in the present calculations range from (the sparsest) mesh with $`20^3`$ nodes to (the thickest) one with $`38^3`$ nodes, i.e. 38 nodes in each direction.

## III Results and discussion

Our results for the energies of the ground state of the two-dimensional $`D^{}`$ center are presented in table 1. This table contains also the corresponding energies for the ground state of the $`D^0`$ center and a comparison with the best results obtained in the literature Larsen1 ; Larsen2 ; Louie . Ref.Larsen1 ; Larsen2  are variational and ref.Louie  represent Monte-Carlo simulations. One can see, that the variational results by Larsen and McCann coincide for weak and intermediate field strengths Larsen2  very well with our results. In contrast to this the results obtained by Louie and Tao PangLouie  differ significantly from these values. The absence of any approximations in our approach and the possibility to reliably evaluate the convergence of our results allow us to conclude, that the results obtained in ref.Louie  overestimate the true values for the binding energy of the ground state of the two-dimensional negative donor considerably. (The reader should note that all digits of the values for our calculated energies given in table 1 are converged i.e. coincide with the exact results). This resolves the discrepancy on the ground state energy of the negative donor present in the literature as demonstrated by the results contained in table 1 particularly for $`\gamma =0`$ but also for nonvanishing field strengths.

For strong magnetic fields $`\gamma >>1`$ table 1 provides also the energies obtained by Larsen and McCann in another variational calculationLarsen1 . Our results both on the total and binding energies show that the ground state of the system is more strongly bound than predicted by ref.Larsen1 . The differences for the energies depend only weakly on the magnetic field strength and are approximately 1.08 a.u. for the total energy and approximately 0.36 a.u. for the binding energy! The significant difference in deviations of our results from those obtained in ref.Larsen1  is due to the fact that the binding energies of the neutral donor in a strong magnetic field are underestimated in ref.Larsen1 . As shown below this leads also to an overestimation of the binding properties of the excited states of the $`D^{}`$ center.

The ground state of the $`D^{}`$ center is the only bound state for $`\gamma =0`$. We have carried out calculations for excited states of both symmetric and antisymmetric character of the spatial wavefunction with respect to the interchange of the coordinates of the electrons (i.e. spin singlet and triplet states). These calculations were implemented both for $`M=0`$ and for $`M0`$. The latter states were investigated also by Larsen and McCannLarsen1 . Without loss of generality we confine ourselves to negative magnetic quantum numbers $`M<0`$. Pairs of states being different only with respect to the sign of $`M`$ possess equal binding energies within the corresponding Landau zones. However, opposite to states with negative $`M`$, states with $`M>0`$ do not belong to the lowest Landau zone. Therefore, it is reasonable to focus on values $`M0`$.

Our calculations show that the ground state of the $`D^{}`$ center is the only spin singlet state which is bound in the presence of magnetic fields. This conclusion coincides with results obtained in ref.Larsen1 . On the otherhand, it follows from our calculations of the triplet states (see table 2) that there are two triplet states that become bound above some corresponding critical values for $`\gamma `$. These are the energetically lowest states for $`M=1`$ and $`M=2`$, respectively. The $`M=1`$ state becomes bound for $`\gamma >0.117`$ (the total energy of both $`D`$ and $`D^{}`$ at $`\gamma =0.117`$ is $`2.057852`$). The state with $`M=2`$ becomes bound for $`\gamma >1.68`$ (the total energy at $`\gamma =1.68`$ is $`2.720788`$).

From table 2 it is evident that the lowest state with magnetic quantum number $`M=3`$ is not bound for any magnetic field strength considered here and probably also for higher field strengths (the binding energy $`E_B=E^DE^D^{}`$ is negative and its absolute value increases with increasing field strength). This result differs from that obtained by Larsen and McCannLarsen1 . They received a positive binding energies for this state for sufficiently strong magnetic fields. Their binding energy for $`\gamma =500`$ is presented in the last row of table 2. The main reason for this discrepancy is the underestimation of the binding energy of the neutral donor in ref.Larsen1  as can be seen from table 1. The data of table 2 together with our results for other excited states (they are unbound) allow us to conclude that the $`D^{}`$ considered possesses only two bound excited states – possessing magnetic quantum numbers $`M=1`$ and $`M=2`$.

Along with the total and binding energies of the two-dimensional negative donor $`D^{}`$ we have calculated some geometrical parameters of its wavefunction that provide additional information about the system. In figure 1 we present the expectation value $`R^2=r_1^2+r_2^2`$ which characterize the extension of the spatial distribution of the electrons in the simplest and most straightforward way as a function of the field strength. First of all we observe that the electronic cloud is shrinking monotonically with increasing field strength for all bound states. This had to expected according to what we know about the behaviour of tightly bound states of few-electron system in strong magnetic fields IvaSchm2001c . The difference of the behavior of $`R^2`$ for the ground ($`M=0`$) and excited ($`M0`$) states for relatively weak magnetic fields is obvious: The ground state is bound for $`\gamma =0`$, its wave function remains localized for all values of $`\gamma `$, and $`R^2`$ possesses a finite limit for $`\gamma 0`$. Furthermore it changes little for $`\gamma 0.1a.u.`$ For the excited states $`R^2`$ is not bounded for $`\gamma 0`$. It changes rapidly with increasing field strength in particular in the weak magnetic field regime. Therefore $`R^2`$ possesses finite values for excited states only due to the presence of the magnetic field. On the otherhand, the dependencies of $`R^2`$ on the magnetic field strength in the high field regime are similar for all considered states since they are dominated by the diamagnetic term $`\gamma (r_1^2+r_2^2)/8`$ of the Hamiltonian. The occurence of a small curvature in the dependence $`R^2(\gamma )`$ for $`M=1`$ between $`\gamma =0.1`$ and $`\gamma =10`$ is due to the influence of internal binding forces of the system, which are not negligibly small compared to the magnetic forces for this range of field strengths.

## IV Brief Summary

We have shown that the two-dimensional negative donor $`D^{}`$ possesses three bound states ie. two bound excited states in the presence of a sufficiently strong magnetic field. The spin singlet ground state is bound for arbitrary field strengths. Our investigation of the ground state clarifies a discrepancy in the literature and confirms previous variational calculations Phelps ; Larsen2 . As a result the energy eigenvalues obtained in ref.Louie  turn out to violate the variational principle and are too low i.e. the corresponding binding energies are too large. For weak and intermediate field strength we obtain good agreement of our total and binding energies compared to those of ref.Larsen2 . In the high field regime, however, a significant lowering of the total energies and raising of the binding energies are obtained within the present investigation. Our particular computational approach allows for an estimate of the difference of our (up to several digits converged) results and the exact ones which consequently allows us to draw definite conclusions on the energies and properties of the donor. A series of calculations for excited states show that two other states become bound with increasing magnetic field strength. They are the lowest excited states with $`M=1`$ ($`\gamma >0.117`$) and $`M=2`$ ($`\gamma >1.68`$). The extension of these states decreases monotonically with increasing field strength.

###### Acknowledgements.

One of the authors (M.V.I.) gratefully acknowledges financial support by the Deutsche Forschungsgemeinschaft. This work was completed during a visit of P.S. to the University of Regensburg whose kind hospitality is appreciated. Financial support by the Graduiertenkolleg ’Nonlinearity and Nonequilibrium in Condensed Matter’ at the University of Regensburg is gratefully acknowledged.


text:
# 1 Introduction

## 1 Introduction

The combined experimental and theoretical study of (ultra-) relativistic heavy-ion collisions has lead to impressive progress in our understanding of the mechanisms underlying the formation and evolution of strongly interacting matter. A key issue concerns the question whether one is indeed producing systems of sufficient macroscopic extent with frequent enough re-interactions among the constituents that justify the use of (thermal-) equilibrium techniques. This, after all, constitutes an essential prerequisite for a meaningful investigation of the phase diagram of QCD in these reactions, including the identification of (the approach towards) phase transitions and Quark-Gluon Plasma (QGP) formation. A successful interpretation of experimental observables has to be based on a comprehensive overall picture of the collision dynamics. For SPS energies, this has been accomplished to a remarkable degree, as evidenced by, e.g., (i) hadro-chemical equilibration as inferred from produced particle ratios, (ii) thermal equilibration with collective transverse expansion, as inferred from transverse momentum spectra combined with two-particle correlation measurements, (iii) a substantial excess of electromagnetic radiation (photons and dileptons) ascribed to thermal emission, (iv) the ’anomalous’ suppression pattern of charmonium states. An important feature of this picture is the distinction between chemical and thermal freezeout, implying a significant duration of an expanding hadronic phase (also essential to explain (ii) and (iii)).

A solid understanding of the majority of observables then allows to go beyond and investigate well-defined deviations. In this talk we will focus on antiproton production, which, upon closer inspection, is not easily understood within the above picture, but, as we argue below, can be reconciled with it.

The paper is organized as follows: in Sect. 2 we briefly recall the essential elements of hadronic equilibration needed to formulate the ’antiproton puzzle’. In Sect. 3 we address the problem within a thermal rate equation, first for the more transparent equilibrium case, then including off-equilibrium effects via its explicit solution. We summarize in Sect. 4.

## 2 Equilibration Times and Antiproton Puzzle

### 2.1 Two Types of Hadronic Freezeout

The notion of chemical and thermal freezeout stages in central collisions of heavy nuclei can be given a clear physical meaning (and applies equally well to hadronic and partonic degrees of freedom). It is based on a hierarchy of time scales related to the underlying cross sections.

Thermal equilibrium in the hadronic fireball can be (locally) maintained as long as elastic collisions are able to keep up with its expansion (dilution) rate. This criterium can be quantified as

$$R(\tau )\underset{\tau }{\overset{\mathrm{}}{}}\frac{d\tau ^{}}{\tau _{th}(\tau ^{})}R_{fo}$$
(1)
with $`R_{fo}1/2`$ (almost independent of collision energy from SIS/BEVALAC to SPS) as extracted from hydrodynamic fits to single-particle spectra . The thermal relaxation time for a hadron species $`i`$ is given by

$$\left(\tau _{th,i}\right)^1=\underset{h}{}\sigma _{ih}v_{rel}\rho _h$$
(2)
in terms of the density $`\rho _h`$ of hadrons which serve as scattering centers, and the pertinent (thermally averaged) cross section $`\sigma _{ih}`$ ($`v_{rel}`$: relative velocity between hadron $`i`$ and $`h`$). Standard (resonance dominated) hadronic reactions, such as $`\pi \pi \rho \pi \pi `$, $`\pi KK^{}\pi K`$, $`\pi N\mathrm{\Delta }\pi N`$ (or even nonresonant $`NN`$ scattering), lead to typical values around $`\sigma 50`$ mb. Even for a rather dilute hadron gas of normal nuclear matter density ($`\rho _h=\rho _0=0.16`$ fm-3) one finds a thermal equilibration time as short as 2 fm/c. A simplified freezeout criterium, which does not explicitly involve the expansion time scale, is given by the condition that the (thermally averaged) mean free path becomes equal to the system size, $`\lambda _{MFP}r_{FB}`$ ($`r_{FB}`$: fireball radius). Both criteria, when applied to central heavy-ion collisions at full SPS energy ($`E_{lab}=158`$ AGeV), indicate thermal freezeout temperatures of about $`T_{th}=110120`$ MeV. This is consistent with combined experimental analysis of (light-) hadron $`p_t`$-spectra and HBT correlations , which can be schematically summarized by (nonrelativistic) slope parameters $`T_{slope,i}=T_{th}+m_iv_t^2`$ with an average transverse flow velocity $`v_t`$ 0.4-0.5 $`c`$ at thermal freezeout.

Chemical equilibrium is maintained by inelastic collisions which change the particle composition of the system. Corresponding hadronic cross sections, e.g., for $`\pi \pi K\overline{K},\rho \rho `$, are typically much smaller than elastic ones, $`\sigma _{inel}1`$ mb after thermal averaging. Consequently, chemical equilibration times are around $`\tau _{chem}100`$ fm/c, well above possible lifetimes of hadronic fireballs created in high-energy nuclear collisions. If anything, the chemical composition of the hadronic phase must frozen early in the evolution. At full SPS energy essentially all hadron abundances can be accommodated by a common chemical freezeout at $`(T_{ch},\mu _B^{ch})(170,260)`$ MeV . This result is, in fact, difficult to understand in purely hadronic scenarios and has been interpreted as ’prehadronic’ flavor equilibration, possibly in a QGP.

A particular consequence of the sequential freezeouts that will be important below is the effective conservation of particle numbers for species that are not subject to strong decays (e.g., $`\pi `$, $`K`$, $`\eta `$; after all, the abundances at thermal freezeout have to agree with the measured ones). In statistical mechanics language this is expressed by finite meson-chemical potentials, which, under SPS conditions, reach appreciable values of $`\mu _\pi ^{th}`$ 60-80 MeV, $`\mu _K^{th}`$ 100-130 MeV (similar for $`\mu _\eta ^{th}`$), etc. Resonances that are regenerated via strong interactions acquire the chemical potentials of their decay products, i.e., $`\mu _\mathrm{\Delta }=\mu _N+\mu _\pi `$, $`\mu _\rho =2\mu _\pi `$, $`\mu _K^{}=\mu _K+\mu _\pi `$, etc.111Strongly decaying resonances with narrow widths ($``$ 50 MeV or so) are in between the two cases. However, for the vector mesons $`\omega `$ and $`\varphi `$, e.g., it turns out that the two limiting cases – full chemical equilibrium ($`\mu _\omega =3\mu _\pi `$ and $`\mu _\varphi =2\mu _k`$) versus chemical off-equilibrium (individual $`\mu _\omega `$ and $`\mu _\varphi `$ adjusted to their conserved number) – give numerically very similar results. We also note that, if reactions of the type $`\omega \pi \pi \pi `$ are relevant, the $`\omega `$ chemical potential would be reduced.

### 2.2 Antiproton Systematics

An enhanced production of antiprotons in central heavy-ion collisions (over the extrapolation from $`pp`$ or $`pA`$ reactions) has been among suggested signatures for QGP formation  (which is expected to facilitate copious antiquark production). The experimental results for $`\overline{p}/p`$ ratios at AGS , SPS  and RHIC  are, however, consistent with the chemical freezeout systematics at the respective energies . In Pb(158 AGeV)+Pb, one finds  $`\overline{p}/p=(5.5\pm 1)`$ %, to be compared to $`\overline{p}/p=\mathrm{exp}(2\mu _N^{ch}/T_{ch})5`$%. Moreover, the $`\overline{p}/p`$ ratio does not exhibit a pronounced centrality dependence , which is also borne out of hadro-chemical model fits .

Despite this apparent agreement it is nevertheless puzzling why the number of antiprotons should to be ’frozen in’ early in the hadronic phase: large $`\overline{p}p`$ annihilation cross section of around 50 mb at the relevant thermal energies ($`E_p^{th}=m_p+\frac{3}{2}kT`$), together with the sizable baryon densities under SPS conditions, seem to imply a chemical freezeout for antibaryons at a much lower temperature (i.e., later time) than $`T_{ch}170`$ MeV; e.g., for $`T=150`$ MeV, where according to thermal fireball calculations  $`\rho _B0.75\rho _0`$, one obtains for the chemical relaxation time,

$$\tau _{\overline{p}}^{chem}=\frac{1}{\sigma _{p\overline{p}}^{ann}(s)v_{rel}\rho _B},$$
(3)
a value of about 3 fm/c, which is significantly smaller than the fireball lifetime. A naive evaluation of the $`\overline{p}/p`$ ratio using thermal freezeout conditions, $`(T_{th},\mu _N^{th})(120,400)`$ MeV , yields $`\mathrm{exp}(2\mu _N^{th}/T_{th})0.1`$ %, a factor $``$ 50 below the measured value! This constitutes the ’antiproton puzzle’.

There have been attempts to resolve this puzzle within sequential scattering (transport) simulations . The latter are, in principle, suited for a realistic description of the late (dilute) hadronic stages in heavy-ion collisions. Nevertheless, extra assumptions, such as a strong in-medium suppression (’shielding’) of the annihilation cross section  or an enhanced production via an increased string tension in the early phases , had to be invoked to avoid or compensate annihilation losses in the hadron gas. However, a notorious problem within these approaches is the treatment of multi-body collisions (3 or more), which are usually neglected, thus violating detailed balance. For the case at hand, these are multi-meson fusion reactions, e.g., $`n_\pi \pi p\overline{p}`$ with an average of $`n_\pi `$ 5-6 pions in the incoming channel. In the following we will show that the inclusion of these reactions, treated within a thermal rate equation, is indeed capable of regenerating a significant number of antiprotons . The notion of detailed balance alone, however, is not enough. A second crucial ingredient  is the over-saturation of pion phase-space, encoded in nonzero pion-chemical potentials, which enhances the equilibrium abundance of antiprotons substantially (see also ref.  for a recent transport approach to the problem).

## 3 Thermal $`\overline{p}`$-Production

### 3.1 Rate Equation and Chemical Potentials

In thermal equilibrium, the (net) rate per unit 4-volume for producing antiprotons via $`n_\pi `$-pion fusion and its back-reaction, $`n_\pi \pi p\overline{p}`$, is given in terms of the corresponding invariant matrix element $`_{n_\pi }`$ as (see, e.g., ref. )

$`_{\overline{p}}^{th}={\displaystyle d^3\stackrel{~}{k}_pd^3\stackrel{~}{k}_{\overline{p}}\underset{i=1}{\overset{n_\pi }{}}d^3\stackrel{~}{k}_i\delta ^{(4)}(K_{tot})|_{n_\pi }|^2\{z_pz_{\overline{p}}\mathrm{e}^{\frac{E_p+E_{\overline{p}}}{T}}z_\pi ^{n_\pi }\mathrm{e}^{\underset{i=1}{\overset{n_\pi }{}}\frac{\omega _i}{T}}\}}.`$ (4)
The thermal factors are in Boltzmann approximation, which allows to factorize the fugacity factors $`z_i=\mathrm{exp}(\mu _i/T)`$. The rate equation can then be written in the compact form

$$_{\overline{p}}^{th}=\frac{d\rho _{\overline{p}}}{dt}\frac{\rho _{\overline{p}}}{\tau _{\overline{p}}^{chem}}\left(\frac{z_\pi ^{n_\pi }}{z_pz_{\overline{p}}}1\right),$$
(5)
which represents a first order differential equation in the antiproton density $`\rho _{\overline{p}}`$. The relaxation time $`\tau _{\overline{p}}^{chem}`$ is given in eq. (3); the use of the total baryon density $`\rho _B`$ therein (rather than just the proton density $`\rho _p`$, which is much smaller) is based on the assumption that the annihilation cross section of antiprotons on other baryons is similar to the one on protons. Eq. (5) furthermore presumes that $`\overline{p}`$ annihilation on excited resonances yields (on average) as many additional pions as the resonance itself decays into (e.g., $`\overline{p}`$ annihilation on a $`\mathrm{\Delta }(1232)`$ gives one pion more than on a nucleon). This leads to a cancellation of additional fugacity factors multiplying $`z_\pi ^{n_\pi }`$ and $`z_p`$ in eq. (5). The use of the total annihilation cross section in $`\tau _{\overline{p}}^{chem}`$ implies a sum over all possible $`n_\pi `$-pion final states ($`s\overline{s}`$ production is suppressed by the OZI rule). Again, this automatically incorporates mesonic resonance states such as $`\rho `$ and $`\omega `$ to the extent that they contribute with a pion-fugacity factor according to their pion final states. In a heavy-ion environment this is obviously satisfied for $`\rho `$ mesons, but to a good approximation also for $`\omega `$ (and even $`\eta `$) mesons, cf. footnote 1.

Under SPS conditions, where the relative amount of antiprotons in the system is small, baryon and pion densities are not significantly altered by the reactions under consideration. The rate equation can thus be considered with the anitproton density as its only variable. In equilibrium one has $`_{\overline{p}}^{th}=0`$, so that the antiproton fugacity can be determined as

$$z_{\overline{p}}^{eq}=z_pz_\pi ^{n_\pi }.$$
(6)
Here, as in eq. (5), the averaging of $`z_\pi ^{n_\pi }`$ is over the initial (or final) pion numbers. Since the pion fugacities enter with large powers, a reasonable estimate of $`z_{\overline{p}}`$ requires a reliable description of the final state multiplicities in $`p\overline{p}`$ annihilation, which we will turn to in the following section.

### 3.2 Properties of $`\overline{p}p`$ Annihilation

Fig. 4 shows the experimental results for pion multiplicities in $`p\overline{p}`$ annihilation at rest (see, e.g., the review  and references therein). The distribution can be reasonably well reproduced by a Gaussian $`P(n_\pi )`$ with an average $`n_\pi =5.06`$ and a width $`\delta n_\pi =0.9`$.

In a thermal heat bath, the annihilation processes occur, on average, at somewhat larger CMS energies, $`\sqrt{s}2E_p^{th}`$. Experimental data on the $`s`$-dependence of $`n_\pi `$ are displayed in Fig. 4. A linear increase accurately describes $`n_\pi (s)`$ in the relevant (thermal) energy range. At the same time, also the width of the Gaussian moderately increases. Put together, the averaged multi-pion fugacity at given temperature takes the form

$$z_\pi ^{n_\pi }(T)=\underset{n_\pi =2}{\overset{n_\pi ^{max}}{}}P(n_\pi )\mathrm{e}^{n_\pi \mu _\pi /T},$$
(7)
where for any practical purposes $`n_\pi ^{max}=9`$.

For the treatment of off-equilibrium effects (cf. Sect. 3.4) one also needs to account for the energy-dependence of the total annihilation cross section, which figures into the equilibration time, eq. (3). Its decrease with increasing $`s`$ can be fitted by the empirical formula $`\sigma _{p\overline{p}}^{ann}(s)=(40p_{lab}^{0.5}+24p_{lab}^{1.1})`$ mb.

### 3.3 $`\overline{p}`$ Thermodynamics at SPS Energies

Inserting the pion-fugacity factor, $`z_\pi ^{n_\pi }(T)`$ from eq. (7) into eq. (6), and using the temperature-dependent pion- and nucleon-chemical potentials as displayed in Fig. 2, results in an equilibrium $`\overline{p}/p`$ ratio for central Pb(158 AGeV)+Pb collisions as represented by the solid lines in Fig. 6.

The enclosed band approximately reflects the uncertainties in determining the pion-chemical potential, associated with the population of higher resonance states (especially those with relatively long lifetime) in the late stages of the hadronic evolution. Irrespective, the strong enhancement over the equilibrium value without the inclusion of finite $`\mu _\pi `$’s (dashed line in Fig. 6) is obvious.

In Fig. 6, the chemical relaxation time is compared to the fireball lifetime, indicating that the antiproton abundance might be driven out of equilibrium before thermal freezeout. This will be quantitatively assessed in the following section.

### 3.4 Off-Equilibrium Evolution

For an explicit solution of the rate equation describing the evolution of the antiproton number in an expanding hadron gas we rely on a hierarchy of time scales $`\tau ^{therm}<\tau _{\overline{p}}^{chem}\tau _\pi ^{chem}`$. This implies an equilibrium $`\overline{p}`$ density

$$\rho _{\overline{p}}^{eq}=z_\pi ^{n_\pi }\frac{d^3k}{(2\pi )^3}\mathrm{exp}[(E_{\overline{p}}\mu _{\overline{p}})/T].$$
(8)
The finite system size in a heavy-ion collision requires the inclusion of volume expansion effects . Using $`dN_{\overline{p}}=_{\overline{p}}^{th}V_{FB}(t)dt`$ where $`N_{\overline{p}}=\rho _{\overline{p}}V_{FB}`$ is the antiproton number and $`V_{FB}(t)`$ the expanding fireball volume , one obtains

$$\frac{d\rho _{\overline{p}}}{dt}=\rho _{\overline{p}}\left(\frac{1}{\tau _{\overline{p}}(t)}+\frac{1}{V_{FB}(t)}\frac{dV_{FB}}{dt}\right)+\rho _{\overline{p}}^{eq}(t)\frac{1}{\tau _{\overline{p}}(t)}\rho _{\overline{p}}(t)L(t)+G(t)$$
(9)
with gain and loss terms $`G(t)`$ and $`L(t)`$, respectively. Eq. (9) is an ordinary first order differential equation which has the solution

$$\rho _{\overline{p}}(t)=\phi (t)\left(\rho _{\overline{p}}(t_0)+\underset{t_0}{\overset{t}{}}𝑑t^{}G(t^{})/\phi (t^{})\right),\phi =\mathrm{exp}[\underset{t_0}{\overset{t}{}}𝑑t^{}L(t^{})].$$
(10)
As shown in Fig. 7, the off-equilibrium results confirm the importance of both the pion-fugacities as well as the backward reaction of multi-pion fusion in maintaining an antiproton abundance close to the observed value. The effective freezeout temperature $`T_{\overline{p}}^{chem}`$ turns out to be close to 130 MeV (slightly higher if one assumes a, e.g., 50% smaller in-medium annihilation cross section).

## 4 Summary

We have shown that antiproton production in central heavy-ion collisions at SPS energies can be reconciled with the observed value corresponding to standard chemical freezeout: annihilation losses in the subsequent hadron gas phase are compensated by the back-reaction of multipion fusion, reinforced by the build-up of large pion-chemical potentials. We have corroborated our earlier findings for the equilibrium case by an explicit solution of the underlying rate equation. An interesting test of the regeneration mechanism could be provided by the recently suggested balance function technique , which has significant sensitivity to the time of (pair) production especially for heavy particles (such as nucleons). The off-equilibrium framework is more involved under RHIC conditions, as $`\rho _B`$ is of the same order as $`\rho _{\overline{B}}`$, so that feedback effects on the baryon densities cannot be neglected. This will be addressed in future work, as it promises valuable insights into hadro-production at collider energies.

ACKNOWLEDGMENT
One of us (RR) would like to thank the organizers for the invitation and a very pleasant and informative meeting. This work has been supported by the US Department of Energy under grant number DE-FG0288ER40388.