text:
# 1 Introduction

## 1 Introduction

The Bose-Einstein correlations (BEC) effect has a quantum-mechanical origin. It arises from the requirement to symmetrise the wave function of a system of two or more identical bosons. It was introduced into particle reactions leading to multi-hadron final states as the GGLP effect  in the study of the $`\pi ^+\pi ^+`$ and $`\pi ^{}\pi ^{}`$ systems. The distributions of the opening angle between the momenta in pairs of like-sign pions were shifted towards smaller values compared to the corresponding distributions for unlike-sign pairs. A related effect was exploited earlier in astronomy  to measure the radii of stars.

In high energy physics, for example $`\mathrm{e}^+\mathrm{e}^{}`$ collisions at LEP, a quantitative understanding of the BEC effect allows tests of the parton fragmentation and hadronisation models. This would in turn help in achieving a more precise measurement of the W boson mass and better knowledge of several Standard Model (SM) observables . The fragmentation models presently used are those of strings and clusters implemented, respectively, in the JETSET  and HERWIG  Monte Carlo generators.

Numerous studies of BEC in pairs of identical bosons already exist, see for example . Due to the experimental difficulties in photon and $`\pi ^0`$reconstruction, only very few studies  exist for BEC in $`\pi ^0`$pairs, even though they offer the advantage of being free of final state Coulomb corrections.

The string model predicts a larger BEC strength or chaoticity and a smaller effective radius of the emitting source for $`\pi ^0`$pairs compared to $`\pi ^\pm `$ pairs while the cluster fragmentation model predicts the same source strength and size . However, neither model of primary hadron production has a mechanism to allow BEC between $`\pi ^0\mathrm{s}`$produced in different strong decays. The string model prediction is a consequence of electric charge conservation in the local area where the string breaks up. Similar expectations can be derived if the probabilities in the string break-up mechanism are interpreted as the squares of quantum mechanical amplitudes . A small difference between $`\pi ^\pm `$ pairs and $`\pi ^0`$pairs is also expected from a pure quantum statistical approach to Bose-Einstein symmetry . In addition, based on isospin invariance, suggestions exist on how to relate BEC in the pion-pair systems i.e. $`\pi ^0\pi ^0,`$ $`\pi ^\pm \pi ^\pm ,`$ and $`\pi ^+\pi ^{}`$ and how to extend it to $`\pi ^\pm \pi ^0`$. The L3 collaboration has recently reported  that the radius of the neutral-pion source may be smaller than that of charged pions, $`R_{\pi ^\pm \pi ^\pm }R_{\pi ^0\pi ^0}=(0.150\pm 0.075(\mathrm{stat}.)\pm 0.068(\mathrm{syst}.))\mathrm{fm},`$ in qualitative agreement with the string fragmentation prediction.

This paper presents a study of BEC in $`\pi ^0`$pairs using the full hadronic event sample collected at centre-of-mass energies at and near the $`\mathrm{Z}^0`$ peak by the OPAL detector at LEP from 1991 to 1995. This corresponds to about four million hadronic $`\mathrm{Z}^0`$ decays. A highly pure sample of $`\pi ^0`$mesons is reconstructed using the lead-glass electromagnetic calorimeter. The correlation function is obtained after accounting for purity and resonant background. It is parametrised with a static picture of a Gaussian emitting source .

## 2 Selection of hadronic $`𝐙^\mathrm{𝟎}`$ decays

A full description of the OPAL detector can be found in . The sub-detectors relevant to the present analysis are the central tracking detector and the electromagnetic calorimeter. The central tracking detector consists of a silicon micro-vertex detector, close to the beam pipe, and three drift chamber devices: the vertex detector, a large jet chamber and surrounding $`z`$-chambers111The OPAL coordinate system is defined so that the $`z`$ axis is in the direction of the electron beam, the $`x`$ axis points towards the centre of the LEP ring, and $`\theta `$ and $`\varphi `$ are the polar and azimuthal angles, defined relative to the $`+z`$\- and $`+x`$-axes, respectively. In cylindrical polar coordinates, the radial coordinate is denoted $`r`$.. In combination, the three drift chambers sitting inside a solenoidal magnetic field of 0.435 T yield a momentum resolution of $`\sigma _{p_t}/p_t\sqrt{0.02^2+(0.0015p_t)^2}`$ for $`|\mathrm{cos}(\theta )|<0.7`$, where $`p_t`$ (in GeV) is the transverse momentum with respect to the beam axis. The electromagnetic calorimeter detects and measures the energies and positions of electrons, positrons and photons for energies above 0.1 GeV. It is a total absorbing calorimeter, and is mounted between the coil and the iron yoke of the magnet. It consists of 11704 lead-glass blocks arranged in three large assemblies (the barrel that surrounds the magnet coil, and two endcaps) which together cover 98% of the solid angle. The intrinsic energy resolution is $`\sigma _E/E5\%/\sqrt{E},`$ where $`E`$ is the electromagnetic energy in GeV.

Standard OPAL selection criteria are applied to tracks and electromagnetic clusters . Tracks are required to have at least 20 measured points in the jet chamber, a measured momentum greater than $`0.1\mathrm{GeV}`$, an impact parameter $`|d_0|`$ in the $`r\varphi `$ plane smaller than 2 cm, a $`z`$ position at the point of closest approach to the origin in the $`r\varphi `$ plane within 25 cm of the interaction point, and a measured polar angle with respect to the beam axis greater than $`20^{}.`$ Electromagnetic clusters are required to have an energy greater than $`0.1\mathrm{GeV}`$ if they are in the barrel part of the detector (i.e. $`\mathrm{cos}\theta 0.82`$) or greater than $`0.3\mathrm{GeV}`$ if they are in the endcap parts. Hadronic $`\mathrm{Z}^0`$ decays are selected by requiring for each event more than 7 measured tracks, a visible energy larger than 60 $`\mathrm{GeV}`$ and an angle larger than $`25^{}`$ and smaller than $`155^{}`$ between the calculated event thrust  axis and the beam axis. The visible energy is the energy sum of all detected tracks, electromagnetic clusters not associated to tracks and electromagnetic clusters associated to tracks after correcting for double counting. A sample of 3.1 million $`\mathrm{Z}^0`$ hadronic decays is selected for which the total background, consisting mainly of $`\tau `$ pairs, is less than $`1\%`$ and is neglected throughout the analysis.

Detector effects and detection efficiencies for the spectra of $`\pi ^0`$pairs are evaluated using eight million Monte Carlo hadronic $`\mathrm{Z}^0`$ decays. Events are generated using the JETSET 7.4 program, tuned to reproduce the global features of hadronic events as measured with the OPAL detector , with the BEC effect explicitly switched off. Samples generated with the HERWIG 5.9 program without the BEC effect are used for comparison. The generated events were passed through a full simulation of the OPAL detector  and were analysed using the same reconstruction and selection programs as were applied to the data.

## 3  Reconstruction of $`\pi ^0`$ mesons

For the selected event sample, neutral pions are reconstructed from photon pairs. Photon reconstruction is performed in the barrel part of the electromagnetic calorimeter where both the photon reconstruction efficiency and the energy resolution are good. The procedure of  which resolves photon candidates in measured electromagnetic clusters is used. It employs a parametrisation of the expected lateral energy distribution of electromagnetic showers. It is optimised to resolve as many photon candidates as possible from the overlapping energy deposits in the electromagnetic calorimeter in a dense environment of hadronic jets. The purity of the photon candidate sample is further increased using a likelihood-type function  that associates to each photon candidate a weight $`w`$ for being a true photon. Photon candidates with higher $`w`$ are more likely to be true photons.

All possible pairs of photon candidates are then considered. Each pair was assigned a probability $`P`$ for both candidates being correctly reconstructed as photons. This probability is simply the product of the $`w`$ -weights associated with the two candidates. The combinatorial background consists of a mixture of three components: (i) wrong pairing of two correctly reconstructed photons, (ii) pairing of two fake photons and (iii) pairing of one correctly reconstructed photon with a fake one. Choosing only photon pairs with high values of $`P`$ leaves combinatorial background mostly from component (i).

The $`\pi ^0`$ reconstruction efficiency and purity are illustrated in Figure 1 for different cuts on $`P.`$ The efficiency is defined as the ratio of the number of correctly reconstructed $`\pi ^0\mathrm{s}`$over the number of generated $`\pi ^0\mathrm{s}`$, and the $`\pi ^0`$purity is defined as the ratio of signal over total entries in a photon-pair mass window between 100 and 170 MeV.

## 4  Selection of $`\pi ^0`$Pairs

The average number of $`\pi ^0\mathrm{s}`$produced in $`\mathrm{Z}^0`$ decays has been measured  to be $`9.76\pm 0.26,`$ which is reproduced by our Monte Carlo simulations. This leads to about $`45`$ possible $`\pi ^0`$pairings per event. Considering only $`\pi ^0`$candidates with $`P>0.1`$ (i.e. 17% efficiency and 36% purity), we reconstruct at the detector level 4.7 $`\pi ^0`$candidates on average per event. This leads to about $`8`$ pairings among which only 1 pair on average is really formed by true $`\pi ^0\mathrm{s}`$. Here, the detector level means that detector response, geometrical acceptance and photon reconstruction efficiency are taken into account. Therefore, the $`\pi ^0`$pair sample is background dominated and the study of $`\pi ^0`$pair correlations or invariant mass spectra is subject to very large background subtraction. Monte Carlo must be used to predict both the shape and amount of background to be subtracted, leading to large systematic errors in the measurements of the BEC source parameters.

To avoid this, the $`\pi ^0`$selection criteria are tightened. We select $`\pi ^0\mathrm{s}`$which have a momentum above 1 GeV. This cut reduces the fraction of fake $`\pi ^0\mathrm{s}`$. In addition, it removes $`\pi ^0\mathrm{s}`$produced by hadronic interactions in the detector material for which the Monte Carlo simulation is not adequate. The probability $`P`$ associated to each $`\pi ^0`$candidate is required to be greater than 0.6. In the case where a photon can be combined in more than one pair, only the pair with the highest probability is considered as a $`\pi ^0`$candidate. Among the events with four or more reconstructed photon candidates, only those leading to a possible $`\pi ^0`$pair with four distinct photon candidates are retained for further analysis. Events with six or more photon candidates leading to more than two $`\pi ^0`$candidates are rejected. They represent about 10% of the retained sample and would increase the sensitivity to unwanted resonance signals if they were not rejected. Figure 2 shows the photon pair mass, $`M_{2\gamma },`$ for the selected events. The average purity of the $`\pi ^0`$sample is 79% in the mass window between 100 and 170 MeV. The background is estimated directly from data by a second order polynomial fit to the side bands of the peak and by Monte Carlo simulation. The two background estimations yield compatible results and the Monte Carlo reproduces correctly the data. The superimposed curves are not the result of a fit to the data, but smoothed histograms of the Monte Carlo expectations for signal and background normalised to the total number of selected hadronic $`\mathrm{Z}^0`$ decays.

A clear $`\pi ^0`$pair signal is obtained as shown in Figure 3 where the two values of $`M_{2\gamma }`$ are shown for the retained events. A $`\pi ^0`$pair is considered as a signal candidate if both values of $`M_{2\gamma }`$ are within the mass window between 100 and 170 MeV. The average $`\pi ^0`$pair signal purity is $`60\%`$ and the Monte Carlo simulation describes the data well. Kinematic fits were made, constraining the mass of pairs of photon candidates to the $`\pi ^0`$mass, with the assumption that the photons come from the primary interaction vertex. Monte Carlo studies showed that this gives a 26% improvement in the resolution of the $`\pi ^0`$momentum.

## 5 The BEC Function

The correlation function is defined as the ratio,

$$C(Q)=\frac{\rho (Q)}{\rho _0(Q)},$$
(1)
where $`Q`$ is a Lorentz-invariant variable expressed in terms of the two $`\pi ^0`$four-momenta $`p_1`$ and $`p_2`$ via $`Q^2=(p_1p_2)^2,`$ $`\rho (Q)=(1/N)\mathrm{d}N/\mathrm{d}Q`$ is the measured $`Q`$ distribution of the two $`\pi ^0\mathrm{s}`$and $`\rho _0(Q)`$ is a reference distribution which should, in principle, contain all the correlations included in $`\rho (Q)`$ except the BEC. For the measurement of $`\rho _0(Q),`$ we consider the two commonly used methods :

* Event Mixing: Mixed $`\pi ^0`$pairs are formed from $`\pi ^0\mathrm{s}`$belonging to different $`\mathrm{Z}^0`$ decay events in the data. To remove the ambiguity on how to mix events, we select two-jet events having a thrust value $`T>0.9,`$ i.e. well defined back-to-back two-jet events. The thrust axes of the two events are required to be in the same direction within ($`\mathrm{\Delta }\mathrm{cos}\theta \times \mathrm{\Delta }\varphi `$)= ($`0.05\times 10^{}`$). Mixing is then performed by swapping a $`\pi ^0`$from one event with a $`\pi ^0`$from another event. To avoid detection efficiency problems arising from different detector regions, swapping of two pions is performed only if they point to the same region of the electromagnetic barrel detector within ($`\mathrm{\Delta }\mathrm{cos}\theta \times \mathrm{\Delta }\varphi `$)= ($`0.05\times 10^{}`$). With this procedure, we start with two hadronic $`\mathrm{Z}^0`$ events each having two $`\pi ^0`$candidates and can end up with between zero and four pairs of mixed $`\pi ^0`$candidates. The $`Q`$ variable is then calculated for each of the mixed pairs. If the contributions from background are removed or suppressed, this method offers the advantage of being independent of Monte Carlo simulations, since $`C(Q)`$ can be obtained from data alone.
* Monte Carlo Reference Sample: The $`\rho _0`$ distribution is constructed from Monte Carlo simulation without BEC. The Monte Carlo is assumed to reproduce correctly all the other correlations present in the data, mainly those corresponding to energy-momentum conservation and those due to known hadron decays. In order to be consistent with the first method, the cut $`T>0.9`$ is also applied for both data and Monte Carlo.

In the following, the distributions $`\rho (Q)`$ and $`\rho _0(Q)`$ are measured from the same sample of selected events. The mixing technique is used as the main analysis method and the Monte Carlo reference technique is applied only for comparison.

## 6 The Measured BEC Function and Background Contribution

The correlation function, $`C(Q),`$ corresponds experimentally to the average number of $`\pi ^0`$pairs, corrected for background, in the data sample divided by the corresponding corrected average number in the reference sample. Thus, we can write

$$C(Q)=\frac{\rho (Q)}{\rho _0(Q)}=\frac{\rho ^m(Q)\rho ^b(Q)}{\rho _0^m(Q)\rho _0^b(Q)},$$
(2)
where $`\rho ^m`$ and $`\rho _0^m`$ are the measured values, and $`\rho ^b`$ and $`\rho _0^b`$ are the corresponding corrections for background contributions. For both the numerator and denominator, the background consists mainly of $`\pi ^0`$pairs in which one or both $`\pi ^0`$ candidates are fake.

The background distributions $`\rho ^b`$ and $`\rho _0^b`$ are obtained from the Monte Carlo information. These background distributions can also be obtained from data using a side band fit to the projected spectra of the two-dimensional $`M_{2\gamma }`$ distributions (see Figure 3) in each 400 MeV interval of the measured $`Q`$ variable. The resulting background distributions are correctly reproduced by Monte Carlo. However, for the smaller $`Q`$ intervals as used in this analysis, i.e. 100 MeV, the side band fit is subject to large statistical fluctuations, so the Monte Carlo distributions have to be used.

In the region of interest where the BEC effect is observed, $`Q<700`$ MeV, pion pairs from particle (resonance) decays could mimic the effect. The relevant decays are: $`\mathrm{K}_\mathrm{s}^0\pi ^0\pi ^0,`$ $`\mathrm{f}_0(980)\pi ^0\pi ^0,`$ and $`\eta \pi ^0\pi ^0\pi ^0`$ with branching ratios of 39%, 33% and 32% respectively. Pion pairs from $`\eta `$ decay contribute only to the region $`Q<315`$ MeV. According to Monte Carlo studies, the number of reconstructed $`\mathrm{K}_\mathrm{s}^0`$ in the 2 $`\pi ^0`$channel is very small. Furthermore, the hypothesis that each $`\pi ^0`$originates from the primary vertex, as used in the kinematic fits (Section 4), does not apply. This is an advantage for this analysis since the $`\mathrm{K}_\mathrm{s}^0`$ peak is flattened, making its effect on the $`Q`$ distribution negligible. The Monte Carlo estimates of this particle decay backgrounds are included in the distribution $`\rho ^b(Q)`$, adjusting the rate of individual hadrons to the LEP average  where necessary.

For our analysis we select $`\pi ^0`$candidates with momentum greater than 1 GeV. This is dictated by the observation of correlations at small $`Q`$ even for Monte Carlo events generated without any BEC effect. Indeed, as shown in Figure 4, a clear BEC-type effect is visible in the correlation function obtained from Monte Carlo events without BEC for different low cuts on $`\pi ^0`$momentum. Using Monte Carlo information, we find that these correlations are mainly caused by $`\pi ^0\mathrm{s}`$originating from secondary interactions with the detector material. They would constitute an irreducible background to the BEC effect if low momentum $`\pi ^0\mathrm{s}`$are considered in the analysis. This effect vanishes for $`\pi ^0`$momenta greater than 1 GeV.

We rely on Monte Carlo simulation only to define the appropriate momentum cut (i.e. 1 GeV) which completely suppresses the effect of soft pions produced in the detector material, rather than relying on its prediction for the exact shape and size of this effect. The reason is that, in contrast to charged pions where the measured track information can be used to suppress products of secondary interactions in the detector material, the neutral pions have to be assumed to originate from the main interaction vertex. Furthermore, with this assumption the kinematic fits (Section 4) bias the energy of soft pions emitted in the detector material towards larger values since the real opening angle between the photons is larger (vertex closer to the calorimeter) than the assumed one.

With the above selection criteria, the composition of the selected $`\pi ^0`$pair sample is studied using Monte Carlo simulations. According to the string fragmentation model implemented in JETSET, the selected sample consists of about 97.9% of mixed pion-pairs from different hadron decays, 2% of pairs belonging to the decay products of the same hadron and only 0.1% prompt pairs from the string break-ups. Similarly, using the cluster fragmentation model implemented in HERWIG, the selected sample consists of 97% of pairs from different hadron decays, 2.3% belonging to the decay products of the same hadron and only 0.7% originating directly from cluster decays. It is worth mentioning that even if the direct pion pairs from string break-up (JETSET) or cluster decays (HERWIG) were all detected and accepted by the analysis procedure, they would be diluted in combination with other pions and would constitute only a marginal fraction ($`<1\%`$ ) of the total number of reconstructed $`\pi ^0`$pairs. Thus, our analysis has no sensitivity to direct pion pairs originating from string break-up or cluster decay.

## 7 Results

The correlation distribution $`C(Q)`$ (Eq. (2)) is parametrised using the Fourier transform of the expression for a static sphere of emitters with a Gaussian density (see e.g. ):

$$C(Q)=N(1+\lambda \mathrm{exp}(R^2Q^2))(1+\delta Q+ϵQ^2).$$
(3)
Here $`\lambda `$ is the chaoticity of the correlation \[which equals zero for a fully coherent (non-chaotic) source and one for a chaotic source\], $`R`$ is the radius of the source, and $`N`$ a normalisation factor. The empirical term, $`(1+\delta Q+ϵQ^2),`$ accounts for the behaviour of the correlation function at high $`Q`$ due to any remaining long-range correlations. The $`C(Q)`$ distribution for data is shown in Figure 5 as the points with corresponding statistical errors, and the smooth curve is the fitted correlation function in the $`Q`$ range between 0 and 2.5 GeV. A clear BEC enhancement is observed in the low $`Q`$ region of the distribution. The parameters are determined to be:

$$\begin{array}{cc}\lambda =& 0.55\pm 0.10,\hfill \\ R=& (0.59\pm 0.08)\mathrm{fm},\hfill \\ N=& 1.10\pm 0.08,\hfill \\ \delta =& (0.14\pm 0.05)\mathrm{GeV}^1,\hfill \\ ϵ=& (0.07\pm 0.03)\mathrm{GeV}^2,\hfill \end{array}$$

where the quoted errors are statistical only and the $`\chi ^2/\mathrm{ndf}`$ of the fit is 14.7/19.

The distribution $`C(Q)`$ obtained for Monte Carlo events generated with no BEC is shown as a histogram in the same figure. It shows that there is no residual correlation at low $`Q`$ and indicates that the observed enhancement is present in the data only. The dashed-line histogram of Figure 5 represents the correlation function obtained from data but before the subtraction, using the Monte Carlo estimates, of pairs from the decay products of the same hadron, indicating that these contributions have only a minor influence on the measured parameters. In addition, the correlation function constructed with background $`\pi ^0`$pairs does not show any enhancement at low $`Q`$(not shown). Here, background $`\pi ^0`$pairs are defined as pairs for which one or both of the $`\pi ^0\mathrm{s}`$are outside the mass window 100-170 MeV, i.e. these are likely to be fake $`\pi ^0`$candidates.

The second method, which uses the MC reference sample, yields the following results:

$$\begin{array}{cc}\lambda & =0.50\pm 0.10,\hfill \\ R& =(0.46\pm 0.08)\mathrm{fm}.\hfill \end{array}$$

These results are quoted for comparison only. We choose to quote the results obtained with the event mixing method since they are much less dependent on details of the Monte Carlo modelling.

The string model predicts a smaller source radius and a larger chaoticity in the BEC effect for $`\pi ^0`$pairs than for $`\pi ^\pm `$ pairs, while the cluster model predicts no difference. These predictions hold only for prompt boson pairs produced directly from the string or cluster decays. According to our Monte Carlo simulations, we have no sensitivity to these pairs.

## 8  Systematic Uncertainties

Potential sources of systematic error are investigated. In each case the effect on the parameters $`R`$ and $`\lambda `$ and their deviations with respect to the standard analysis are estimated. The results are summarised in Table 1.

* Bin width resolution: After the kinematic fits (Section 4), the resolution on the invariant mass of two pions, or on the variable $`Q,`$ is approximately 60 MeV. We have chosen a bin width of 100 MeV for the fit to the measured $`C(Q)`$ distribution. This bin width is varied from 100 MeV to 80 MeV and to 120 MeV.
* Fit range: The low end of the fit range is set to start at $`Q=350`$ MeV (fourth bin) The high end of the fit range is changed to stop at $`Q=2`$ GeV.
* Effect of hadron decays:  To estimate the effect of the $`\pi ^0`$pairs from the same resonance decay on the measured BEC parameters, the estimated contribution is varied by $`\pm 10\%`$ which represents the typical error on the measured individual hadron rates . In order to investigate the dependence of the measured parameters $`R`$ and $`\lambda `$ on the $`\pi ^0`$momentum cut, the analysis is repeated for $`\pi ^0`$momenta larger than 1.2 GeV.
* Analysis procedure:  The analysis is repeated for several variations of the selection criteria.

+ 1) The $`\pi ^0`$selection mass window is changed from 100-170 MeV to 110-165 MeV ( increases the $`\pi ^0`$purity by 5%).
    + 2) The probability for $`\pi ^0`$selection is changed from 0.6 to 0.5 (reduces the $`\pi ^0`$purity by 5%).
    + 3) The thrust value for two-jet events is changed from 0.9 to 0.85 and to 0.92 (changes the overall event sample size by $`\pm 5\%`$).
    + 4) The factor $`1+\delta Q+ϵQ^2`$ is replaced by $`1+\delta Q.`$
    + 5) $`\pi ^0`$from different events are mixed if they point to the same region of the detector within $`(\mathrm{\Delta }\mathrm{cos}\theta \times \mathrm{\Delta }\varphi )=(0.10\times 15^{})`$ instead of $`(0.05\times 10^{}).`$

The final systematic errors are obtained from quadratically adding the deviations from the central value. Thus,

$$\begin{array}{cc}\lambda & =0.55\pm 0.10\pm 0.10,\hfill \\ R& =(0.59\pm 0.08\pm 0.05)\mathrm{fm}.\hfill \end{array}$$

## 9 Conclusions

We have observed Bose-Einstein correlations of $`\pi ^0`$pairs produced in hadronic $`\mathrm{Z}^0`$ decays. Assuming a Gaussian shape for the source, we obtain $`\lambda =0.55\pm 0.10\pm 0.10`$ for the chaoticity parameter and $`R=(0.59\pm 0.08\pm 0.05)\mathrm{fm}`$ for the radius. In order to construct a reference sample with the event mixing method, this analysis is restricted to well defined back-to-back two-jet events. Furthermore, in order to remove $`\pi ^0\mathrm{s}`$not originating from the primary interaction vertex the considered momentum phase space is restricted to $`p_{\pi ^0}>1`$ GeV. The measured value of the source radius is smaller than our former value , $`R=(1.002\pm 0.016_{0.096}^{+0.023})\mathrm{fm},`$ obtained for charged pions for which the measured track parameters allowed access to lower momenta and where the reference sample was constructed with unlike-sign pion pairs. However, the value is compatible with the LEP inclusive average , $`R=(0.74\pm 0.01\pm 0.14)\mathrm{fm}`$, for charged pions. Pions from strong decays constitute the dominant part of our sample of reconstructed $`\pi ^0`$pairs. We have no sensitivity to test the string or cluster model predictions concerning differences between neutral and charged pion pairs. We deduce that Bose-Einstein correlations exist between $`\pi ^0`$ pairs in which each $`\pi ^0`$is a strong decay product of a different hadron.

## Acknowledgements

We particularly wish to thank the SL Division for the efficient operation of the LEP accelerator at all energies and for their close cooperation with our experimental group. In addition to the support staff at our own institutions we are pleased to acknowledge the
Department of Energy, USA,
National Science Foundation, USA,
Particle Physics and Astronomy Research Council, UK,
Natural Sciences and Engineering Research Council, Canada,
Israel Science Foundation, administered by the Israel Academy of Science and Humanities,
Benoziyo Center for High Energy Physics,
Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) and a grant under the MEXT International Science Research Program,
Japanese Society for the Promotion of Science (JSPS),
German Israeli Bi-national Science Foundation (GIF),
Bundesministerium für Bildung und Forschung, Germany,
National Research Council of Canada,
Hungarian Foundation for Scientific Research, OTKA T-029328, and T-038240,
The NWO/NATO Fund for Scientific Reasearch, the Netherlands.


format: text

text:
# Suppression of Rabi oscillations for moving atoms

## I Introduction

The coupling between the center-of-mass motion of the atom and localized laser fields is a central topic in quantum optics, and leads to many important effects and applications in cooling, trapping, deflection, and isotope separation experiments. The objective of this paper is to point out one further consequence of such a coupling: the suppression of Rabi oscillations and the related possibility of controlled projection onto internal pure states by the motion of the atom from a laser-free region into a laser-illuminated region.

For a general two-level system with an interaction-picture Hamiltonian of the form

$$H=\frac{\mathrm{}}{2}\left(\genfrac{}{}{0pt}{}{0}{\mathrm{\Omega }}\genfrac{}{}{0pt}{}{\mathrm{\Omega }}{2\delta }\right),$$
(1)
a state $`𝚿\left(\genfrac{}{}{0pt}{}{\mathrm{\Psi }^{(1)}}{\mathrm{\Psi }^{(2)}}\right)`$, beginning in the ground state $`|1=\left(\genfrac{}{}{0pt}{}{1}{0}\right)`$ at $`t=0`$, will evolve according to

$`\mathrm{\Psi }^{(1)}`$ $`=`$ $`e^{i\delta t/2}\left[\mathrm{cos}\left({\displaystyle \frac{\mathrm{\Omega }^{}t}{2}}\right)i{\displaystyle \frac{\delta }{\mathrm{\Omega }^{}}}\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Omega }^{}t}{2}}\right)\right],`$
$`\mathrm{\Psi }^{(2)}`$ $`=`$ $`{\displaystyle \frac{ie^{i\delta t/2}\mathrm{\Omega }}{\mathrm{\Omega }^{}}}\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Omega }^{}t}{2}}\right),`$ (2)
where

$$\mathrm{\Omega }^{}=(\mathrm{\Omega }^2+\delta ^2)^{1/2}$$
(3)
and $`\mathrm{\Omega }`$ is the usual Rabi frequency.

There is currently much interest in the dynamics of two-level systems governed by the Hamiltonian of Eq. (1) plus an additional time-dependent driving term, in particular to determine the Rabi oscillation suppression for critical values of the external field parameters. In the so called dynamical localization effect, for example, the system remains in one of the two levels Sacchetti01 . Another Rabi-oscillation-suppression effect has been described recently for the scattering of atoms by a standing light wave EFYS02 . In this case, the partial contributions from different diffraction angles lead to different oscillation periods and thus to inhomogeneous broadening.

In this paper we describe a different type of oscillation suppression due to the motion of the system rather than to an additional time dependent field, or to a diffraction effect. Explicitly, we shall consider the motion of a two-level atom with transition frequency $`\omega `$ in one dimension, with a classical electric field (laser of frequency $`\omega _L`$) illuminating the half axis $`x>0`$ perpendicularly. Without damping, and in an interaction picture for the internal degrees of freedom, the Hamiltonian can be written in the form

$$H=\widehat{p}^2/2m+\frac{\mathrm{}}{2}\left(\genfrac{}{}{0pt}{}{0}{\mathrm{\Omega }\mathrm{\Theta }(\widehat{x})}\genfrac{}{}{0pt}{}{\mathrm{\Omega }\mathrm{\Theta }(\widehat{x})}{2\delta }\right),$$
(4)
where $`\mathrm{\Theta }`$ is the step function, $`\mathrm{\Omega }`$ plays the role of a laser-atom coupling constant, $`\delta =\omega _L\omega `$ is the detuning, and hats denote operators whenever confusion is possible with a corresponding $`c`$-number. In particular, $`\widehat{p}`$ is the momentum operator for the $`x`$ direction. If one takes damping due to photon emissions into account, it can be shown, by means of the quantum jump approach Hegerfeldt91 , that in three dimensions the atomic time development between emissions may be described by an (effective, non-hermitian) “conditional” Hamiltonian $`H_\mathrm{c}`$,

$`H_\mathrm{c}`$ $`=`$ $`\widehat{𝐩}^2/2m+{\displaystyle \frac{\mathrm{}}{2}}\mathrm{\Omega }\mathrm{\Theta }(\widehat{x})\{|21|e^{ik_L\widehat{y}}+\mathrm{h}.\mathrm{c}.\}`$ (5)
$``$ $`{\displaystyle \frac{\mathrm{}}{2}}(2\delta +i\gamma )|22|,`$
where $`\gamma `$ is the Einstein coefficient of level 2, i.e. its decay rate or inverse life time, and $`\widehat{𝐩}`$ is the momentum operator in three dimensions (3D). The factor $`e^{ik_L\widehat{y}}`$ takes into account the spatial dependence of the laser coupling. A Hamiltonian of the form of Eq. (4) is obtained from Eq. (5) by neglecting spontaneous emissions, i.e. setting $`\gamma =0`$, and assuming in addition that the atomic wave packet is centered at $`y=0`$ and satisfies $`k_L\mathrm{\Delta }y1`$, at least for some time, so that the exponentials can be dropped and a one dimensional kinetic term suffices. This approximation and its limitations will be further commented on Appendix A.

The conditional Hamiltonian $`H_\mathrm{c}`$ is closely related to waiting times between photon emissions Cohen-Dalibard85 . Indeed, let an atomic state $`|\mathrm{\Psi }(0)`$ be prepared at time $`0`$. Then it can be shown by means of the quantum jump approach Hegerfeldt91  that

$$P_0(t)e^{iH_\mathrm{c}(t)/\mathrm{}}|\mathrm{\Psi }(0)^2$$
(6)
gives the probability of no emission until time $`t`$ and $`\mathrm{\Pi }(t)P_0^{}(t)`$ is the probability density for the first photon. After an emission the atom has to be reset, here to its ground state, and then the conditional time development resumes. In this way one can simulate an emission sequence and the corresponding Bloch (master) equation for the atom.

As we will show, for moving atoms there are Rabi oscillations both in time as well as in space. In the temporal case one asks for the probability, $`P_2(t)`$, of finding the atom in the excited state, without regard to its spatial position. For the simple model in Eq. (4) without damping, this can be easily calculated and Figure 1 depicts a suppression for this case. The solid line shows the population of the excited state $`P_2`$ versus $`t`$ for $`\delta =0`$ and for a wave packet which is prepared at $`t=0`$ with negligible negative momentum components as a minimum-uncertainty-product Gaussian, with the atom in the ground state, and far from the laser. At sufficiently large times, and for the parameters we have chosen, the wave packet is completely within the laser region, but there are no Rabi oscillations. Instead, $`P_2`$ increases monotonically, and saturates at $`0.5`$; this might look surprising at first sight: there is no classical averaging (the state is quantum and coherent), there is no dispersion related broadening, there is no time dependent term in the Hamiltonian, there is not even decay that could smooth or suppress the oscillation.

Basically, for measurements in time the effect has two different regimes and explanations depending on the incident energy: at very low kinetic energy, the suppression is due to the formation of a pure non-oscillating internal state of the atom in the laser region, whereas at intermediate energies it may be understood in terms of a simple semiclassical approximation. Yet at higher energies the effect disappears, and the expected Rabi oscillations may be seen more and more clearly for increasing velocities, as illustrated by the dashed and dotted lines of Figure 1.

Spatial Rabi oscillations can occur in the emission probability at different atomic positions. If $`P_2(x,t)`$ denotes the probability density to find the atom in its excited state at the position $`x`$ at time $`t`$ and if one defines $`I(x)`$ by

$$I(x)𝑑tP_2(x,t),$$
(7)
then $`\gamma I(x)dx`$ gives the number of emitted photons (per atom) with the atom at $`dx`$, and thus $`I(x)`$ is proportional to the photon intensity for the atom at position $`x`$ foot1 . Neglecting damping, the expression on the rhs can be calculated with Eq. (4), and an example is given in Fig. 2, where the spatial Rabi oscillation is completely absent in one of the two cases depicted. In contrast to the temporal case, spatial oscillations are only suppressed by the quantum mechanism, i.e., at very low kinetic energies.

In the next four sections we shall characterize, neglecting damping due to spontaneous emissions, several aspects of the Rabi oscillation suppression due to atomic motion: the basic theory (Section II), the effect of the velocity regimes in the Rabi oscillations for measurements in time (Section III), the effect of detuning and the velocity in the internal atomic states (Section IV), and the suppression of spatial Rabi oscillations (Section V). The effect of damping will be taken into account in Section VI.

## II Basic theory

### II.1 Exact results

Here we shall present the theory required to describe the interaction of the moving atom with the perpendicular laser. We first neglect damping and solve the eigenvalue equation for the Hamiltonian of Eq. (4) subject to the condition that the atom impinges on the laser beam from the left in the ground state. The eigenvalues are $`E=\mathrm{}^2k^2/2m`$ and the eigenfunctions are denoted by $`𝚽_k\left(\genfrac{}{}{0pt}{}{\varphi ^{(1)}}{\varphi ^{(2)}}\right)`$, with $`k>0`$.

For $`x0`$, $`𝚽_k`$ is of the form

$$𝚽_k(x)=\frac{1}{\sqrt{2\pi }}\left(\genfrac{}{}{0pt}{}{e^{ikx}+R_1e^{ikx}}{R_2e^{iqx}}\right),x0,k>0,$$
(8)
where

$$q^2=k^2+2m\delta /\mathrm{}$$
(9)
and $`\mathrm{Im}(q)0`$. For large enough (negative) detuning, $`k<|2m\delta /\mathrm{}|^{1/2}`$, $`q`$ becomes purely imaginary and the reflected wave for the excited state decays exponentially.

In the laser region, let $`|\lambda _+`$ and $`|\lambda _{}`$ be the eigenstates of the matrix $`\frac{1}{2}\left(\genfrac{}{}{0pt}{}{0}{\mathrm{\Omega }}\genfrac{}{}{0pt}{}{\mathrm{\Omega }}{2\delta }\right)`$ corresponding to the eigenvalues $`\lambda _\pm `$. One easily finds

$`\lambda _\pm `$ $`=`$ $`{\displaystyle \frac{1}{2}}[\delta \pm \mathrm{\Omega }^{}],`$ (10)
$`|\lambda _\pm `$ $`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{1}{\frac{(\delta \pm \mathrm{\Omega }^{})}{\mathrm{\Omega }}}}\right).`$ (11)
Note that $`|\lambda _\pm `$ have not been normalized.

For $`x0`$, one can write $`𝚽_k`$ as a superposition of $`|\lambda _\pm `$,

$$\sqrt{2\pi }𝚽_k(x)=C_+|\lambda _+e^{ik_+x}+C_{}|\lambda _{}e^{ik_{}x}.$$
(12)
\[The mathematical solutions $`e^{ik_\pm x}`$ are not included because they correspond to negative momenta or increasing exponentials.\]

From the eigenvalue equation $`H𝚽_k=E_k𝚽_k`$ it follows that

$$k_\pm ^2=k^22m\lambda _\pm /\mathrm{}=k^2+\frac{m}{\mathrm{}}(\delta \pm \mathrm{\Omega }^{}),$$

with Im $`k_\pm 0`$, and from the continuity of $`𝚽_k(x)`$ at $`x=0`$ one obtains, with $`|1=\left(\genfrac{}{}{0pt}{}{1}{0}\right)`$ and $`|2=\left(\genfrac{}{}{0pt}{}{0}{1}\right)`$,

$`1+R_1`$ $`=`$ $`C_+1|\lambda _++C_{}1|\lambda _{}`$
$`R_2`$ $`=`$ $`C_+2|\lambda _++C_{}2|\lambda _{}.`$
Similar relations result from the continuity of $`𝚽_k^{}(x)`$ at $`x=0`$, yielding

$`C_+`$ $`=`$ $`2k(q+k_{})\lambda _{}/D,`$
$`C_{}`$ $`=`$ $`2k(q+k_+)\lambda _+/D,`$
$`R_1`$ $`=`$ $`[\lambda _+(q+k_+)(kk_{})\lambda _{}(q+k_{})(kk_+)]/D,`$
$`R_2`$ $`=`$ $`k(k_{}k_+)\mathrm{\Omega }/D,`$ (13)
where

$$D=(k+k_{})(q+k_+)\lambda _+(k+k_+)(q+k_{})\lambda _{}.$$
(14)
Thus Eq. (12) becomes, in components and for $`x0`$,

$`\varphi _k^{(1)}(x)`$ $`=`$ $`{\displaystyle \frac{2k}{\sqrt{2\pi }D}}\{(q+k_{})\lambda _{}e^{ik_+x}(q+k_+)\lambda _+e^{ik_{}x}\},`$
$`\varphi _k^{(2)}(x)`$ $`=`$ $`{\displaystyle \frac{k\mathrm{\Omega }}{\sqrt{2\pi }D}}\{(q+k_{})e^{ik_+x}(q+k_+)e^{ik_{}x}\}.`$ (15)
Time dependent wave packets incident on the laser region with positive momentum components and with the atom in the ground state may be formed by linear superposition DEHM02 ,

$$𝚿(x,t)=_0^{\mathrm{}}𝑑ke^{i\mathrm{}k^2t/2m}𝚽_k(x)\stackrel{~}{\psi }(k),$$
(16)
where $`\stackrel{~}{\psi }(k)`$ is the wave number amplitude at time $`t=0`$ corresponding to the freely moving packet.

The conditional Hamiltonian $`H_c`$ for the one-dimensional case with damping can be treated in a similar way, and this is done in Appendix B.

### II.2 Different approximation regimes

Depending on the incident velocity different approximations are applicable. For large velocities a semiclassical approximation is sufficient whereas for slow atoms it is essential to retain the quantum nature of the translational motion.

#### II.2.1 Fast atoms: semiclassical approximation

For “high” kinetic energy, $`k^2|m(\delta \pm \mathrm{\Omega }^{})/\mathrm{}|`$, we may approximate the wave numbers in the wave function by

$`k_\pm `$ $``$ $`k+{\displaystyle \frac{m}{2\mathrm{}k}}[\delta \pm \mathrm{\Omega }^{}],`$ (17)
$`q`$ $``$ $`k+{\displaystyle \frac{m\delta }{k\mathrm{}}},`$ (18)
and therefore,

$`\varphi _k^{(1)}(x)`$ $``$ $`{\displaystyle \frac{e^{ikx}}{\sqrt{2\pi }}}\{e^{i\delta t/2}[\mathrm{cos}\left({\displaystyle \frac{m\mathrm{\Omega }^{}x}{2\mathrm{}k}}\right)`$
$``$ $`i{\displaystyle \frac{\delta }{\mathrm{\Omega }^{}}}\mathrm{sin}\left({\displaystyle \frac{m\mathrm{\Omega }^{}x}{2\mathrm{}k}}\right)]\},`$
$`\varphi _k^{(2)}(x)`$ $``$ $`{\displaystyle \frac{e^{ikx}}{\sqrt{2\pi }}}\left\{e^{i\delta t/2}\left[{\displaystyle \frac{i\mathrm{\Omega }}{\mathrm{\Omega }^{}}}\mathrm{sin}\left({\displaystyle \frac{m\mathrm{\Omega }^{}x}{2\mathrm{}k}}\right)\right]\right\}.`$ (19)
\[The relation with the Raman-Nath approximation is discussed in Appendix C.\] Up to the normalization factor this reveals a “spatial Rabi oscillation” that may be associated with the temporal one by the change of variable $`t=xm/k\mathrm{}`$, the time that a classical particle with momentum $`k\mathrm{}`$ needs to travel from the origin to $`x`$. With this change the terms in curly brackets correspond exactly to the internal state amplitudes of the atom at rest, see Eq. (2), while the multiplying plane wave denotes an undisturbed center of mass motion with fixed momentum.

Moreover in this regime reflection is negligible,

$$R_1\left(\frac{m\mathrm{\Omega }}{4\mathrm{}k^2}\right)^2,R_2\frac{m\mathrm{\Omega }}{4\mathrm{}k^2}.$$
(20)
A further simplification is to assume that the interferences among different momenta are either negligible or average out. Then we may quite accurately approximate the center of mass motion of a wave packet with high kinetic energy by a purely classical free motion, whereas the internal state is treated quantum mechanically. More specifically, in this approximation the wave packet width and spreading are taken into account by using an ensemble of semiclassical atoms with their center of mass distributed in classical phase space according to the initial quantum Wigner function $`W(q_0,p_0;t=0)`$. \[For free motion, the quantum and “classically evolved” Wigner functions coincide at all times if initially identical.\] To each of these atoms we associate an internal state, that starts to oscillate once the center of mass crosses $`x=0`$ at $`q_0m/p_0`$. The expectation values in the semiclassical approximation are then computed by averaging over the ensemble of semiclassical atoms (translationally classical, internally quantum). In particular, the probability to find the atom in the excited state at a time $`t`$ is given by

$`P_2(t)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑q_0{\displaystyle _0^{\mathrm{}}}𝑑p_0W(q_0,p_0;t=0)`$
$`\times `$ $`\left({\displaystyle \frac{\mathrm{\Omega }}{\mathrm{\Omega }^{}}}\right)^2\mathrm{sin}^2[\mathrm{\Omega }^{}(t+q_0m/p_0)/2]\mathrm{\Theta }(t+q_0m/p_0),`$
where we have assumed that $`t>0`$, that the wave packet at $`t=0`$ is far from the laser, and that all atoms have positive momentum.

#### II.2.2 Slow atoms

At the critical wavenumber $`k=k_c`$,

$$k_c[m(\mathrm{\Omega }^{}\delta )/\mathrm{}]^{1/2},$$
(22)
the de Broglie wavelength for the atomic motion $`\lambda _{dB}(k)=2\pi /k`$ becomes equal to the “spatial period” of the Rabi oscillation in Eq. (19), or “Rabi wavelength” $`\lambda _R(k)=2\pi k/k_c^2`$. For smaller $`k`$ the quantum aspect of the atomic motion cannot be ignored, and $`k_{}`$ becomes purely imaginary,

$$k_{}ik_c,$$
(23)
whereas $`k_+`$ remains real. As a consequence, $`e^{ik_{}x}`$ evanesces, and the contribution of $`|\lambda _{}`$ beyond $`k_c^1`$ (the penetration length for this component) vanishes,

$`\varphi _k^{(1)}(x)`$ $``$ $`{\displaystyle \frac{C_+}{\sqrt{2\pi }}}e^{ik_+x},x>k_c^1`$ (24)
$`\varphi _k^{(2)}(x)`$ $``$ $`{\displaystyle \frac{C_+(\delta +\mathrm{\Omega }^{})}{\sqrt{2\pi }\mathrm{\Omega }}}e^{ik_+x},x>k_c^1.`$
This results in a pure, and non-oscillating internal state, as discussed later in Section IV.

## III Temporal Rabi oscillations: velocity effects

In addition to the separation between “classical” and “quantum” regimes at $`k_c`$, see Eq. (22), the classical regime may be also subdivided into intermediate and high velocities, depending on whether or not the temporal Rabi oscillations are suppressed. (We shall see later than for spatial oscillations this later subdivision does not apply.)

### III.1 High velocities: $`k>k_R`$

Suppose that the center of mass motion of the wave packet can be reproduced with an ensemble of classical atoms, as discussed in the previous section, and that they enter into the laser region in a very short time compared to the the Rabi period $`T_R2\pi /\mathrm{\Omega }^{}`$ (sudden entrance). In that case they will oscillate in phase and the Rabi oscillations will be seen, as in the dotted line of Figure 1. For Gaussian packets, the transition (average) wave number between adiabatic and sudden entrance regimes may be identified by equating $`T_R`$ and a measure of the duration of the wave packet passage across the origin, $`t_p5\mathrm{\Delta }_xm/k\mathrm{}`$, where $`\mathrm{\Delta }_x`$ is the wave packet width (square root of variance) at the peak’s passage across $`x=0`$, and the factor 5 is rather arbitrary but is only intended to give an estimate,

$$k_R\frac{5\mathrm{\Delta }_xm\mathrm{\Omega }^{}}{2\pi \mathrm{}}.$$
(25)
### III.2 Intermediate velocities: $`k_c<k<k_R`$

In the intuitive language suggested by the classical approximation for the center of mass motion, the “mechanism” that explains the suppression of temporal Rabi oscillations at intermediate velocities is the averaging of the Rabi oscillations of the individual atoms forming the ensemble when the entrance time for the whole ensemble is long compared to the Rabi period. Since each atom enters at a different instant in the laser region, the phase of the oscillation will be also different. This is the case shown in Fig. 3, where the peaks of the densities of the ground and excited state alternate with a spatial period given by $`\mathrm{}kT_R/m`$.

Figure 4 shows the “visibility” $`V_t`$ obtained from maxima (max) and minima (min) of the probability of the excited state, $`P_2(t)`$, corresponding to large times beyond the initial transient,

$$V_t=\frac{P_2(max)P_2(min)}{P_2(max)+P_2(min)},$$
(26)
versus the incident average velocity for five different values of $`\mathrm{\Omega }`$ and a fixed wave packet width, as well as the estimate given by Eq. (25); again $`\delta =0`$. The visibilities calculated exactly or with the classical approximation are indistinguishable.

### III.3 Low velocity: $`k<k_c`$

Below $`k_c`$ there are no Rabi oscillations either but, according to the discussion of section II.2.2, the reason is not an effective “averaging” (the semiclassical approximation is not valid), but the fact that the atomic internal state formed by the laser does not oscillate at all.

## IV Internal states: Effects of the velocity and detuning

The motion of the atom from the laser-free region to the laser-illuminated region has very different effects on the internal state depending on the incident velocity and the detuning. We shall discuss now the effect of the various velocity regimes in the degree of mixing of the normalized, internal, reduced density operator for the atoms in the laser region, which is defined by its matrix elements

$$\rho _{ij}^L(t)=\frac{_0^{\mathrm{}}𝑑x[\psi ^{(i)}(x,t)]^{}\psi ^{(j)}(x,t)}{_0^{\mathrm{}}𝑑x|\psi ^{(1)}(x,t)|^2+|\psi ^{(2)}(x,t)|^2},$$
(27)
distinguishing between the cases of zero and non-zero detuning.

### IV.1 Zero detuning

(i) For $`k<k_c`$ the populations of both ground and excited states are equal and the overlap in space of the respective wave functions is maximum. After a transient time, necessary for Eq. (24) to apply, $`\rho ^L`$ is finally given by the pure state $`2^{1/2}(|1|2)`$,

$$\rho ^L=\left(\genfrac{}{}{0pt}{}{1/2}{1/2}\genfrac{}{}{0pt}{}{1/2}{1/2}\right),$$
(28)
which is nothing but $`|\lambda _+\lambda _+|`$, now normalized.

(ii) For $`k_c<k<k_R`$ excited and ground components do not overlap, see Fig. 3, so $`\rho ^L`$ tends to a diagonal matrix. The two populations are still equal (there is no global Rabi oscillation because of the semiclassical averaging effect) and $`\rho ^L`$ is 1/2 times the unit matrix, which amounts to the maximum degree of mixing (incoherence) allowed for the reduced internal state. The sharp transition from a pure state to a maximally mixed state around $`k=k_c`$ may be seen in Figure 5.

There is also quite a dramatic effect of $`k_c`$ on the reflection probabilities, which tend to vanish for higher values of $`k`$, see Figures 6 and 7 for $`|R_1|^2`$ and $`|R_2|^2`$ respectively.

(iii) Finally, for $`kk_R`$ there is no classical averaging effect, and the two populations oscillate periodically and coherently in time between 0 and 1. The reduced internal state in the laser region becomes pure again, as Figure 5 shows. Note that the transition from intermediate to high kinetic energies at $`k_R`$ is rather smooth, in particular in comparison with the sharp $`k_c`$-transition.

If the initial atom is prepared in the excited state instead of the ground state, different scattering eigenfunctions have to be used. The calculation follows similar lines and provides different expressions for reflection amplitudes and coefficients $`C_\pm `$. However, the eigenstates $`|\lambda _\pm `$ remain the same, and it is still true that below $`k=k_c`$ only $`k_+`$ is real, so that the same pure state is obtained at low kinetic energy. This means that the low-energy motion of the atom into a laser illuminated region provides a very simple physical mechanism to project any initial internal state, pure or mixed, onto the pure state $`2^{1/2}(|1|2)`$ if $`\delta =0`$. The laser illuminated region may in fact be moved and the atom be at rest with the same effect.

### IV.2 Non-zero detuning

Detuning may be used as a control knob to obtain states with a smaller fraction of excited state. The (unnormalized) pure state obtained at low kinetic energy in the laser region is given for arbitrary $`\delta `$ by

$$𝚿=\left(\genfrac{}{}{0pt}{}{1}{\frac{\delta +\mathrm{\Omega }^{}}{\mathrm{\Omega }}}\right).$$
(29)
As seen in Figure 5, for the intermediate regime $`k_c<k<k_R`$, the classical averaging is not so effective for non zero detuning, so that the state is less mixed. Note also that negative detuning leads to less reflection at very low kinetic energies than positive detuning, see Figs. 6 and 7, and the existence of an additional critical point associated with the transition between imaginary and real values of $`q`$.

## V Suppression of spatial Rabi oscillations

The quantum suppression of Rabi oscillations below $`k_c`$ may also be seen in several space dependent quantities. They are, however, not sensitive at all to the semiclassical mechanism: these quantities oscillate both below and above $`k_R`$ defined in Eq (25).

In Figure 8 we show

$$V_x=\frac{|\varphi _k^{(2)}(x)|^2(max)|\varphi _k^{(2)}(x)|^2(min)}{|\varphi _k^{(2)}(x)|^2(max)+|\varphi _k^{(2)}(x)|^2(min)}$$
(30)
for maxima and minima with respect to $`x`$ evaluated beyond the penetration length of the $`\lambda _{}`$-component. One may see a rather abrupt jump from 0 to 1 at $`k=k_c`$.

A similar jump may be observed for $`I(x)`$ of Eq. (7). For wave packets with components below $`k_c`$, so that $`k_+k_c,k_{}ik_c`$, and for $`x`$ beyond the transient, small-$`x`$ region, the $`t`$ integral can be done to obtain $`I(x)v/v_c^2`$. On the other hand, for $`k`$ above $`k_c`$, $`I(x)`$ becomes

$$I(x)\frac{m}{\mathrm{}}_0^{\mathrm{}}\frac{dk}{k}\mathrm{sin}^2\left(\frac{\mathrm{\Omega }x}{2v}\right)|k|\psi ^{(1)}(0)|^2.$$
(31)
These two limiting cases are depicted in Figure 2.

## VI Inclusion of Damping

In this section we shall see that the suppression of Rabi oscillations can be detected by means of the fluorescence signal of the spontaneously emitted photons. Two cases will be considered: a detection measurement of the first fluorescence photon, which is described by the conditional Hamiltonian, and a measurement where all successive photons are also taken into account, which requires a master equation (optical Bloch equations). In both cases the fluorescence signal is proportional to the population of the excited state $`P_2`$ which can be calculated by means of $`H_c`$ or by the Bloch equations, respectively.

The suppression of Rabi oscillations discussed here should be distinguished from their fading away due to damping and approach to a stationary internal state. The latter is quite well-known for atoms at rest WM  (cf. Figure 9). It is possible to operationally distinguish fading and suppression by a combination of $`\mathrm{\Omega }`$ and $`\gamma `$ that, for the atom at rest, allows to observe several oscillations before reaching the asymptotic stationary population. This corresponds to strong driving conditions, $`\mathrm{\Omega }\gamma `$.

### VI.1 Rabi oscillations in the first-photon distribution

For the one-dimensional case the conditional Hamiltonian $`H_c`$ can be written as

$$H_c=\frac{\widehat{p}^2}{2m}+\frac{\mathrm{}}{2}\left(\genfrac{}{}{0pt}{}{0}{\mathrm{\Omega }\mathrm{\Theta }(\widehat{x})e^{ik_Ly}}\genfrac{}{}{0pt}{}{\mathrm{\Omega }\mathrm{\Theta }(\widehat{x})e^{ik_Ly}}{2\delta i\gamma }\right).$$
(32)
The generalization of Eqs. (11-14) for $`\gamma 0`$ is provided in Appendix B.

For high velocities the entrance of the wave packet in the laser region is sudden, so that $`P_2`$ versus time takes the same form as for the atom at rest from the entrance time on. At lower velocities the oscillation pattern disappears; $`k_R`$ still marks the transition, which is now smoother than when $`\gamma =0`$. Figure 10 shows $`P_2`$ for a velocity below $`k_R`$. The Rabi oscillation suppression is evident, compare with the dashed line of Figure 9.

The transition at $`k_c`$ is illustrated in Fig. 11, which shows the degree of mixing versus the average velocity foot2 . Note the smoothing effect of a non-zero $`\gamma `$.

An observation of the suppression of Rabi oscillations using the first photon could be achieved for a “Lambda” configuration of three atomic levels such that the laser couples two of them, while a third one acts as a sink for the excited state EFYS02  (only one photon may be emitted per atom). The ideal experiment would require a preparation of the atoms, sent one by one, according to a single wave function, and a detector capable of responding to a single photon. The experiment would be repeated many times until a profile similar to the one in Fig. 10 is obtained for the density of photon detection times.

### VI.2 Temporal Rabi oscillations in the photon intensity

For a simple two-level atomic configuration, without a sink level for the excited state, the atom may emit many photons. In this case recoil effects may become important, especially at low incident energies, and lead to spatial broadening. The master equation that describes the density operator $`\rho ^{3D}`$ of the atomic system in three spatial dimensions is given, following Hegerfeldt93 , by

$$\dot{\rho }^{3D}=\frac{i}{\mathrm{}}\{H_c\rho ^{3D}\rho ^{3D}H_c^{}\}+(\rho ^{3D}),$$
(33)
where $`H_c`$ is the conditional 3D Hamiltonian, Eq. (5), and the “reset” or “jump” term $`(\rho ^{3D})`$ takes into account the atomic recoil along the directions weighted by the dipole emission distribution. In momentum representation,

$`𝐩|(\rho ^{3D})|𝐩^{}`$ $`=`$ $`\gamma |11|`$
$`\times `$ $`{\displaystyle 𝑑𝜿P(𝜿)𝐩+\mathrm{}k_L𝜿|\rho _{22}^{3D}|𝐩^{}+\mathrm{}k_L𝜿},`$
where the integral is over all possible photon directions represented by the unit vectors $`𝜿=𝐤/|𝐤|`$, and

$$P(𝜿)=\frac{3}{8\pi }\left(1\frac{𝜿𝐝}{|𝐝|^2}\right)=\frac{3}{8\pi }\mathrm{sin}^2\theta ,$$
(35)
$`\theta `$ being the angle between $`𝜿`$ and the dipole moment $`𝐝`$, which is taken along the $`z`$ direction.

Even though the fluorescence photons and the associated atom recoil may go in any direction, one can take the trace over momentum components $`p_y`$ and $`p_z`$ in the reset term, and reduce the integral to a one dimensional one, see e.g. Hensinger01 , to obtain for the initial atomic-motion direction $`x`$ the expression

$`_x(\rho )`$ $``$ $`\mathrm{Tr}_{y,z}(\rho ^{3D})=|11|\gamma {\displaystyle \frac{3}{8}}`$
$`\times `$ $`{\displaystyle _1^1}𝑑u(1+u^2)p_x+\mathrm{}k_Lu|\rho _{22}|p_x^{}+\mathrm{}k_Lu,`$
where the reduced atomic density operator for the $`x`$-direction is

$$\rho _{ij}\mathrm{Tr}_{y,z}i|\rho ^{3D}|j,i=1,2,j=1,2.$$
(37)
Note that it contains both internal and center of mass information of the atom, so it is different from $`\rho ^L`$ in Eq. (27). Even though, taking the trace over $`y`$ and $`z`$ in the other terms of Eq. (33) does not lead to a closed equation for the reduced one-dimensional density matrix of Eq. (37), a one-dimensional approximation is valid,

$$\dot{\rho }=\frac{i}{\mathrm{}}\{H_c\rho \rho H_c^{}\}+_x(\rho ),$$
(38)
with $`H_c`$ given by Eq. (32), provided that the effect of the two broadening mechanisms in the $`y`$ direction, the standard quantum mechanical wave packet spreading and the atomic recoil, remain small. The time scales of both effects are estimated in Appendix A.

The fluorescence signal is proportional to $`P_2`$, which is now the total population of excited atoms, regardless of whether or not they have emitted. The one dimensional master equation may be solved in principle directly, or, as done here, using the quantum jump technique, i.e., averaging over many “trajectories” Hegerfeldt91 . For each of them the photon detections occur at random instants. Figure 12 shows $`P_2`$ averaged over 10000 trajectories, compared with the solid line of Figure 9, properly shifted, which gives the time evolution of $`P_2`$ for a sudden entrance of a very fast wave packet. The suppression of Rabi oscillations is evident.

The ideal fluorescence measurement would be performed by sending one atom at a time, as described before. In a less demanding experiment, one could prepare a cloud of many atoms of sufficiently low density, and send it towards the laser, measuring the fluorescence signal versus time. The cloud ensemble would lead to suppression of Rabi oscillations at sufficiently low velocities. In this case the incident translational state would not be pure but a true mixture.

The transition at $`k_c`$ would be most easily noticed by the reflection dominant at lower kinetic energies. For the Cs transition we are considering in the numerical examples $`v_c\mathrm{}k_c/m=0.28`$ m/s, which is well above the recoil velocity limit $`\mathrm{}k_L/m=0.35`$ cm/s. The possibility to control the pure state formed in the laser region will depend on the ability to focus the atomic beam along the $`y`$\- direction in the scale of the laser wavelength, due to the $`y`$-dependence of the pure state obtained, see Eq. (48).

### VI.3 Spatial Rabi oscillations

With damping included, $`P_2(x,t)`$ becomes $`x|\rho _{22}(t)|x`$ so that $`I(x)`$ in Eq. (7) becomes

$$I(x)=_{\mathrm{}}^{\mathrm{}}𝑑tx|\rho _{22}(t)|x,$$
(39)
and $`\gamma I(x)`$ is the mean number of photons per atom per unit length.

In a similar way, one may consider the probability density, denoted by $`\gamma I_0(x)`$, for the atomic position when the first photon is emitted. Similar to Eq. (7), $`I_0(x)`$ is given by

$$I_0(x)=_{\mathrm{}}^{\mathrm{}}𝑑t|\psi ^{(2)}(x,t)|^2,$$
(40)
where the time development is given by the conditional Hamiltonian in Eq. (32), and $`\psi ^{(2)}(x,t)`$ denotes the excited-state component. Integrating Eq. (40) over $`x`$ gives the fraction of atoms which emit a photon, so that $`𝑑xI_0(x)`$ tends to one for $`kk_c`$.

Above $`k_c`$, and using the language of the classical approximation, the time when the atom arrives at the laser and starts the Rabi cycle is unimportant for an observation of the spatial dependence of the photon intensity; the only relevant information is the position where the emission takes place. In other words, the suppression will not be visible at the intermediate velocities below $`k_R`$; for an example see Figure 13, where both $`\gamma I_0(x)`$ and $`\gamma I(x)`$ are depicted for the same initial conditions, with no suppression of the oscillations. In the opposite “quantum” case, $`k<k_c`$, the suppression effect is visible in $`I_0(x)`$; for an example see the dashed line of Figure 13. Note the clear distinction between the small-$`x`$ region, where there is still a contribution from $`|\lambda _{}`$, and the larger-$`x`$ region, which depends only on $`|\lambda _+`$. A similar result holds for $`I(x)`$ but, because of the continuation of pumping-emission cycles, there is no exponential decay; see Figure 14 for an example. We may estimate $`I(x)`$ after the transient peak by approximating $`x|\rho _{22}(t)|x`$ in Eq. (39) by its value at $`\gamma =0`$, as in Section V. This gives, with $`kk_c`$, $`I(x)v/v_c^2`$. The exponential decay of $`I_0(x)`$ may also be approximated in the strong driving limit and $`kk_c`$ by retaining the dominant exponentially decaying terms in Eq. (24), $`I_0(x)ve^{\gamma (m/\mathrm{\Omega }\mathrm{})^{1/2}/2}/v_c^2`$.

## VII Discussion

The possibility to reach ultra cold temperatures for atomic motion leads to a number of new interesting quantum phenomena that would remain hidden otherwise. As is well known, the internal level populations of fast atoms exhibit Rabi oscillations when traversing a sufficiently strong traveling laser field. In this paper, however, we have shown that for cold atoms these Rabi oscillations may be suppressed and in fact can be totally absent. In the context of moving two-level atoms, two types of Rabi oscillations can occur, temporal and spatial ones. In the temporal case one observes the population of the upper level in time and the associated photon emission, regardless of the atomic position. In the spatial case one observes where the emissions occur in space, without regard to time.

We have distinguished two regimes for the suppression of temporal Rabi oscillations at low and intermediate kinetic energies and have characterized the transitions quantitatively. Spatial oscillations are only suppressed at low kinetic energies.

Quite generally the effects of slow atomic motion or laser displacement on a two-level qubit are of interest in the design on quantum computers SH03 . Applications of the present effect may be based on the fact that at low kinetic energy the suppression of Rabi oscillations is associated with the projection onto specific internal pure states, which can be controlled via detuning. Purification procedures are essential in quantum data processing CEM99 . Our results are also relevant for the analysis of time scales defining the particle’s traversal through a region of space tqm . They suggest that the excited state population cannot in general be taken as a measure of the time elapsed in the laser region (based solely on a comparison with the Rabi oscillation of the atom at rest) unless the particle’s kinetic energy is sufficiently high.

Even though we have centered the discussion on an atom-laser system, the formalism is equally valid for other two-level systems, and in particular for a localized magnetic field acting on a spin-$`1/2`$ particle, i. e., to the so called “Larmor-clock”, which is usually applied with an additional potential barrier and in the weak coupling limit tqm .

A simplification of our treatment has been the semi-infinite laser-illuminated region. This is not a major obstacle to observe the effect since, due to the low velocities implied, actual finite laser profiles would be enough for the occurrence of many Rabi oscillations of a classically moving atom. For example, a Rabi period for our Cs transition takes 0.036 $`\mu `$s when $`\mathrm{\Omega }=5\gamma `$, whereas an atom with speed ten times $`v_c`$ would only move 0.1 $`\mu `$m in that period.

###### Acknowledgements.

We are grateful to R. F. Snider for very useful comments. This work has been supported by Ministerio de Ciencia y Tecnología (BFM2000-0816-C03-03), UPV-EHU (00039.310-13507/2001), and the Basque Government (PI-1999-28).

## Appendix A Limits of the 1D model

A simple random, classical model will enable us to estimate a maximum time of validity of the one dimensional approximation used in the text. At variance with other random walk approaches dealing with stimulated emission and absorption processes in standing waves ABS81 , the present model focuses on the momentum kicks due to spontaneous photon emissions.

Suppose that the atom suffers random kicks at $`j=1,2,3,\mathrm{},n`$ instants separated by the interval $`\mathrm{\Delta }t`$ and that at $`t=0`$ the atom is at rest, $`y=v=0`$, where $`v`$ is the velocity in the y-direction. If $`\mathrm{\Delta }v_j`$ and $`\mathrm{\Delta }y_j`$ are velocity and position increments in the $`y`$-direction in step $`j`$ and $`v_y`$ and $`y_j`$ the corresponding velocity and position, we have

$$y_n=\underset{j}{}\mathrm{\Delta }y_j=\mathrm{\Delta }t\underset{j}{}v_j$$
(41)
and

$`{\displaystyle \underset{j}{}}v_j`$ $`=`$ $`(\mathrm{\Delta }v_1)+(\mathrm{\Delta }v_1+\mathrm{\Delta }v_2)+\mathrm{}+(\mathrm{\Delta }v_1+\mathrm{}+\mathrm{\Delta }v_n)`$
$`=`$ $`n\mathrm{\Delta }v_1+(n1)\mathrm{\Delta }v_2+(n2)\mathrm{\Delta }v_3+\mathrm{}.+\mathrm{\Delta }v_n,`$
so that

$`y_n^2`$ $`=`$ $`(\mathrm{\Delta }t)^2[n^2(\mathrm{\Delta }v_1)^2+(n1)^2(\mathrm{\Delta }v_2)^2+\mathrm{}+(\mathrm{\Delta }v_n)^2]`$
$`+`$ $`\mathrm{crossed}\mathrm{terms}\mathrm{involving}\mathrm{\Delta }v_i\mathrm{\Delta }v_l,il.`$
When taking the average, the contribution from the crossed terms vanishes,

$`y_n^2`$ $`=`$ $`(\mathrm{\Delta }t)^2[n^2+(n1)^2+\mathrm{}+1](\mathrm{\Delta }v)^2`$
$`=`$ $`{\displaystyle \frac{n(n+1)(2n+1)}{6}}(\mathrm{\Delta }t)^2(\mathrm{\Delta }v)^2,`$
where $`(\mathrm{\Delta }v)^2`$ is the variance of the velocity increments in the $`y`$ direction at any of the $`j`$ instants. For large $`n`$,

$$y^2\frac{n^3}{3}(\mathrm{\Delta }t)^2(\mathrm{\Delta }v)^2.$$
(42)
For the dipole radiation distribution of Eq. (35), and considering that the modulus of each velocity increment is due to the recoil velocity $`\mathrm{\Delta }v=k_L\mathrm{}/m`$, we find $`(\mathrm{\Delta }v_y)^2=2\mathrm{\Delta }v^2/5`$. A critical $`n`$ may be obtained from (42) by imposing, say, $`k_Ly^2^{1/2}=1/10`$. For the transition of Cs atoms we are considering, $`\gamma =33.3\times 10^6`$ s<sup>-1</sup>, $`k_L^1=852/2\pi `$ nm, and for strong driving $`\mathrm{\Delta }t2/\gamma `$; this gives $`n30`$, or approximately 2 $`\mu s`$ of time from the first photon emission.

The other phenomenon to take into account is quantum mechanical dispersion of the wave packet. This may be estimated by assuming $`k_Ly^2^{1/2}=1/10`$ in the formula

$$(\mathrm{\Delta }y)^2=(\mathrm{\Delta }y)_0^2\left[1+\frac{\mathrm{}^2t^2}{4m^2(\mathrm{\Delta }y)_0^4}\right].$$
(43)
For $`k_L(\mathrm{\Delta }y)_0=1/20`$, this gives 3 $`\mu `$s.

## Appendix B Theory with decay

We provide generalizations of Eqs. (11-14) for $`\gamma 0`$ and an arbitrary value of $`y`$ corresponding to the Hamiltonian of Eq. (32):

$$\lambda _\pm =\frac{1}{2}\left[\delta +i\frac{\gamma }{2}\pm \left(\delta ^2\frac{\gamma ^2}{4}+i\delta \gamma +\mathrm{\Omega }^2\right)^{1/2}\right],$$
(44)
$$|\lambda _\pm =\left(\genfrac{}{}{0pt}{}{1}{\frac{2\lambda _\pm }{\mathrm{\Omega }e^{ik_Ly}}}\right),$$
(45)
$$k_\pm ^2=k^2\frac{2m\lambda _\pm }{\mathrm{}},q^2=k^2+\frac{2m}{\mathrm{}}(\delta +i\gamma /2),$$
(46)
where $`k_\pm `$ and $`q`$ have positive imaginary parts. In terms of these variables, the expressions for $`C_\pm `$, $`R_1`$, and $`D`$ are the same as in Eqs. (13) and (14), whereas

$$R_2=k(k_+k_{})\mathrm{\Omega }e^{ik_Ly}/D,$$
(47)
In this one-dimensional approximation the parametric dependence on $`y`$ does not affect the reflection probabilities, $`\lambda _\pm `$, $`q`$, $`k_\pm `$, or $`C_\pm `$, but it gives a phase factor to $`R_2`$, and to $`2|\lambda _\pm `$. A consequence of the later is the formation of different pure states in the laser region at low kinetic energy. In particular, for $`\gamma =0`$,

$$|\lambda _+=\left(\genfrac{}{}{0pt}{}{1}{\frac{\delta +\mathrm{\Omega }^{}}{\mathrm{\Omega }e^{ik_Ly}}}\right).$$
(48)
## Appendix C Relation with the Raman-Nath approximation

It is interesting to compare the approximation of the text ($`x`$ as a variable, $`y`$ as a parameter), with the Raman-Nath (RN) approximation where $`y`$ is the variable, and $`x`$ a parameter, expressed in terms of time through an assumed (classical) linear relation

$$t=mx/k\mathrm{},$$
(49)
where $`k`$ is also a parameter. Consider the atom evolving with the Hamiltonian

$$H_{RN}=\frac{\mathrm{}}{2}\left(\genfrac{}{}{0pt}{}{0}{\mathrm{\Omega }e^{ik_Ly}}\genfrac{}{}{0pt}{}{\mathrm{\Omega }e^{ik_Ly}}{2\delta }\right).$$
(50)
Note the absence of a kinetic term, a basic feature of the RN approximation. This Hamiltonian is easily diagonalized. The eigenvalues are given by Eq. (10) and the eigenvectors take the form of Eq. (45) so that the wave vector for a state initially in the ground state with amplitude $`\psi _0(y)`$ is given by

$`\psi _{RN}^{(1)}(y)`$ $`=`$ $`\psi _0(y)\left[{\displaystyle \frac{\mathrm{\Omega }^{}\delta }{2\mathrm{\Omega }^{}}}e^{i\lambda _+t}+{\displaystyle \frac{\mathrm{\Omega }^{}+\delta }{2\mathrm{\Omega }^{}}}e^{i\lambda _{}t}\right],`$
$`\psi _{RN}^{(2)}(y)`$ $`=`$ $`\psi _0(y)\left[{\displaystyle \frac{\mathrm{\Omega }e^{ik_Ly}}{2\mathrm{\Omega }^{}}}e^{i\lambda _+t}+{\displaystyle \frac{\mathrm{\Omega }e^{ik_Ly}}{2\mathrm{\Omega }^{}}}e^{i\lambda _{}t}\right],`$
or $`𝚿_{RN}=\psi _0𝚽_{RN}`$. Making the substitution of Eq. (49) in $`𝚽_{RN}`$ one finds the same result as $`𝚽_𝐤`$ of Eq. (12) when using the “high kinetic energy” approximations of Eqs. (17,18), except for the plane wave factor $`e^{ikx}/(2\pi )^{1/2}`$,

$$𝚽_k\frac{e^{ikx}}{(2\pi )^{1/2}}𝚽_{RN}.$$
(51)
Physically, the plane wave accounts for the undisturbed translational motion along the incident direction whereas the internal one is described by the RN wave vector $`𝚽_{RN}`$.

In the RN approximation the probability to find the atom in the excited state is given by $`\mathrm{\Omega }^2\mathrm{sin}^2(t\mathrm{\Omega }^{}/2)/\mathrm{\Omega }^2`$. In spite of the neglect of the spatial displacement along $`y`$, there is a possibility of momentum exchange due to the absorption of a laser photon TRC84 . Taking the Fourier transform we find the following normalized momentum distributions for ground and excited state

$`𝒫_1(p_y)`$ $`=`$ $`|p_y|\psi _0|^2,`$ (52)
$`𝒫_2(p_y)`$ $`=`$ $`|p_yk_L\mathrm{}|\psi _0|^2.`$ (53)
Note that the RN approximation cannot describe the Rabi-oscillation-suppression effect at low kinetic energies or the transition region around $`k_c`$ since the longitudinal motion is treated classically. The suppression effect at intermediate energies could be described if an ensemble of “RN atoms” mimicking a wave packet distributed along the longitudinal direction is considered.


text:
# Was the Universe Reionized by Massive Metal Free Stars?

## 1 Introduction

The first stars in the universe (Pop-III) formed out of metal-free gas, relic from the big bang. The bottom-up hierarchy of structure formation in Cold Dark Matter (CDM) models implies that stars started to form in CDM halos just above the cosmological Jeans mass of $`10^5M_{}`$ at redshifts as high as $`z2030`$ (see review by Barkana & Loeb 2001 and references therein). Simulations of metal-free star formation indicate that the first stars had high masses ($`M100M_{}`$) since gas cooling by molecular hydrogen (H2) cannot lower the gas temperature below $`200`$K (Bromm, Coppi, & Larson 2002; Abel, Bryan & Norman 2002). Massive, metal-free stars shine at their Eddington luminosity, $`L_\mathrm{E}M`$, and have roughly constant effective (surface) temperatures of $`10^5`$K and lifetimes of $`3\times 10^6`$yr, independent of their mass (Bromm, Kudritzki, & Loeb 2001). This implies that the number of ionizing photons produced per baryon incorporated into these stars ($`4.5\times 10^4`$) was independent of their mass function, and larger by more than an order of magnitude than the yield of the observed Pop-II (low-mass, metal-rich) stars.

Did the first generation of metal-free stars produce a sufficient total number of ionizing photons to reionize the universe? This question is difficult to answer theoretically. The UV radiation from the massive stars can easily dissociate the fragile intergalactic H2 molecules which are responsible for star formation in mini-halos with virial temperatures $`T_{\mathrm{vir}}10^{2.4}`$K (Haiman, Rees, & Loeb 1997; Haiman, Abel & Rees 2000; Ricotti et al. 2002). Moreover they could evaporate the gas out of the shallow potential wells of the first CDM mini-halos (Barkana & Loeb 1999; Nishi & Tashiro 2000). Finally the hydrodynamic energy released by a single supernova explosion is much larger than the gravitational binding energy of the gas in these halos (Barkana & Loeb 2001). If feedback quenches star formation in mini-halos, then Pop-III star formation will be dominated by more massive halos where H2 can form in dense gas shielded from the UV background (Oh & Haiman 2002; Mackey, Bromm & Hernquist 2002).

Analysis of the first year data from the WMAP satellite suggests a large optical depth to electron scattering $`\tau _{\mathrm{es}}=0.17\pm 0.04`$, implying significant reionization at redshifts as high as $`17\pm 5`$ (Kogut et al. 2003; Spergel et al. 2003). While this result is not surprising in view of earlier detailed calculations of the reionization history by metal-free stars (Wyithe & Loeb 2003; Paper I), it offers an intriguing empirical path for answering the aforementioned question. In this Letter, we argue that metal-free stars are required by the WMAP measurement of $`\tau _{es}`$.

We assume the most recent cosmological parameters obtained through fits to WMAP data (Spergel et al. 2003). These include density parameters of $`\mathrm{\Omega }_m=0.27`$ in matter, $`\mathrm{\Omega }_b=0.044`$ in baryons, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.73`$ in a cosmological constant, a Hubble constant of $`H_0=71\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, an rms amplitude of $`\sigma _8=0.84`$ for mass density fluctuations in a sphere of radius $`8h^1`$Mpc, and a primordial power-spectrum with a power-law index $`n=1`$. Any tilt in the power spectrum at the level that is marginally indicated by the combination of the WMAP data with other data sets (Spergel et al. 2003) would reduce the small-scale power, delay reionization, and strengthen our conclusion.

## 2 Optical Depth to Electron Scattering

Different reionization histories result in different large-scale polarization anisotropies of the microwave background (Zaldarriaga 1997; Kaplinghat et al. 2002). A simple measure of the reionization history is the integrated optical depth to electron scattering $`\tau _{\mathrm{es}}=_0^{10^3}𝑑z|cdt/dz|\sigma _Tn_\mathrm{e}(z)`$, where $`\sigma _T`$ is the Thomson cross-section and $`(dt/dz)^1=H_0(1+z)\sqrt{\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }}`$. The electron number density, $`n_\mathrm{e}`$, depends on the mass filling factors $`Q^\mathrm{m}`$ of the different ionized species in the gas,

$$n_\mathrm{e}=(1+z)^3\left[Q_{\mathrm{H}^+}^\mathrm{m}n_\mathrm{H}^0+Q_{\mathrm{He}^+}^\mathrm{m}n_{\mathrm{He}}^0+2Q_{\mathrm{He}^{++}}^\mathrm{m}n_{\mathrm{He}}^0\right],$$
(1)
where $`n^0`$ denotes the total present-day number densities of hydrogen (H), and helium (He).

## 3 Reionization Histories

We now consider reionization histories under different model assumptions. The histories follow the reionization of hydrogen as well as the single and double reionization of helium in an inhomogeneous IGM (Miralda-Escude et al. 2000). We include helium in our calculations since the harder spectrum of massive Pop-III stars can reionize both hydrogen and helium at high redshift (Oh et al. 2001; Venkatesan, Tumlinson & Shull 2003).

We calculate the emissivity of quasars described by the luminosity function of Wyithe & Loeb (2002a) as well as the emissivity of high redshift zero-metallicity stars and Pop-II stars (the latter is discussed in § 3.1). Our model has two free parameters: (i) the transition redshift, $`z_{\mathrm{tran}}`$, at which Pop-II stars start to dominate the cosmic radiation field of ionizing photons; and (ii) the product of the star formation efficiency in massive galaxies $`f_{}`$ (=stellar mass/total initial baryonic mass) and the escape fraction $`f_{\mathrm{esc}}`$ of ionizing photons from galaxies. In Paper I we presented contours of $`\tau _{\mathrm{es}}`$ as a function of $`z_{\mathrm{tran}}`$ and $`f_{}f_{\mathrm{esc}}`$. We have found that if the star formation efficiency in early low-mass galaxies resembled that in local dwarf galaxies with the same circular velocity, then values of $`f_{}f_{\mathrm{esc}}0.050.2`$ and $`z_{\mathrm{tran}}15`$ lead to $`\tau _{\mathrm{es}}`$ within the range measured by WMAP. Below we summarize the basic ingredients of our model; more details can be found in Paper I.

We treat the IGM as a two-phase medium. In regions of reionized IGM, star formation can only be initiated in halos that collapse above a virial temperature threshold of $`T_{\mathrm{min}}=2.5\times 10^5`$K, corresponding to the cosmological Jeans mass of photo-ionized IGM at $`10^4`$K (e.g., Thoul & Weinberg 1996; Kitayama & Ikeuchi 2000). In the neutral IGM, star formation is limited to halos with virial temperatures above the minimum for efficient gas cooling. If only atomic hydrogen is present, then $`T_{\mathrm{min}}=10^4`$K. However if molecular hydrogen is present then $`T_{\mathrm{min}}=10^{2.4}`$K (Haiman, Thoul, & Loeb 1996; Haiman, Abel, & Rees 2000) and star formation can proceed in much smaller halos. We begin by discussing the case with $`T_{\mathrm{min}}=10^4`$K.

### 3.1 Ionizing Photons from Stars

To compute the comoving densities of stellar ionizing photons ($`n_\gamma `$) we integrate over model spectra for stellar populations to find the total number of ionizing photons ($`N_\gamma `$) emitted per baryon incorporated into stars. We find

$$\frac{dn_\gamma }{dz}=N_\gamma f_{\mathrm{esc}}\frac{dF_{}}{dz}n_\mathrm{b},$$
(2)
where $`n_\mathrm{b}`$ is the number density of baryons, and $`F_{}(z)`$ is the fraction of baryons in the universe at redshift $`z`$ that collapsed, cooled, and formed stars.

A recent analysis of SDSS data (Kauffmann et al. 2002) has shown that in a large sample of local galaxies, the ratio between the total stellar ($`M_{}`$) and dark matter halo ($`M`$) masses scales as $`M^{2/3}`$ for $`M_{}<3\times 10^{10}M_{}`$, but is constant for larger stellar masses. Thus the star formation efficiency within dwarf galaxies is proportional to $`M^{2/3}`$. Since star formation is thought to be regulated by supernova feedback (Dekel & Silk 1986; Dekel & Woo 2002), and the supernova energy output is proportional to $`M_{}`$ while the gravitational binding energy scales as $`Mv_{\mathrm{cir}}^2M^{5/3}`$, more massive galaxies retain a larger fraction of their gas. Using the stellar mass Tully-Fisher relation of Bell & De Jong (2001), we find the threshold circular velocity $`v_{}=176\mathrm{km}\mathrm{s}^1`$ that corresponds to a stellar mass of $`3\times 10^{10}M_{}`$ at $`z=0`$. The mass of a dark matter halo with circular velocity $`v_{}`$ is $`M_{\mathrm{halo}}^{}5\times 10^{10}\left[(1+z)/10\right]^{3/2}M_{}`$. We define $`f_{}`$ as the star formation efficiency in galaxies with circular velocities larger than $`v_{}`$; below $`v_{}`$ the star formation efficiency is $`\eta _{}=f_{}\left(M/M_{\mathrm{halo}}^{}\right)^{2/3}`$. Given $`T_{\mathrm{min}}`$ we find

$$\frac{dF_{}(T_{\mathrm{min}})}{dz}=\frac{d}{dz}\frac{1}{\rho _\mathrm{m}}_{M(T_{\mathrm{min}})}^{\mathrm{}}𝑑Mϵf_{}M\frac{dn_{\mathrm{st}}}{dM},$$
(3)
where $`ϵ=1`$ for $`M>M_{\mathrm{halo}}^{}`$ and $`ϵ=\left(M/M_{\mathrm{halo}}^{}\right)^{2/3}`$ for $`M<M_{\mathrm{halo}}^{}`$, $`\rho _\mathrm{m}`$ is the mean comoving mass density of matter, and $`n_{\mathrm{st}}(>M)`$ is the comoving number density of halos with masses $`>M`$ (Sheth & Thormen 1999). Finally, the value of $`dF_{}/dz`$ in the two phase IGM required for calculation of $`dn_\gamma /dz`$ may then be computed from

$$\frac{dF_{}}{dz}=Q_{\mathrm{H}^+}^\mathrm{m}\frac{dF_{}(2.5\times 10^5\text{K})}{dz}+(1Q_{\mathrm{H}^+}^\mathrm{m})\frac{dF_{}(10^4\text{K})}{dz}.$$
(4)
### 3.2 The Need for Metal-Free stars

We begin with calculations of reionization histories in the absence of metal-free stars at all redshifts. We have calculated reionization histories for several values of $`f_{}f_{\mathrm{esc}}`$ in terms of the volume filling factor $`Q`$ of ionized regions in the low overdensity ($`20`$) IGM (Miralda-Escude et al. 2000; Paper I), and computed the resulting values of $`\tau _{\mathrm{es}}`$.

Recently, Theuns et al. (2002) have argued that if HI reionization ended before $`z=9`$, and the universe remained reionized after that, then adiabatic cooling would have lead to an IGM temperature at $`z=3.5`$, much lower than inferred from the Ly$`\alpha `$ forest. As a result, the most recent overlap of intergalactic HII regions must have occurred at $`z9`$, unless there was significant photo-heating at $`z3.5`$ due to the double reionization of helium. In the absence of metal free stars we find that this constraint requires $`f_{}f_{\mathrm{esc}}0.1`$, resulting in a value of $`\tau _{\mathrm{es}}0.09`$ that is not consistent with the WMAP data. On the other hand, if we neglect the constraint that HII overlap occurred at $`z9`$, we find that in the absence of metal-free stars, reionization histories resulting in $`\tau _{\mathrm{es}}0.13`$ imply HII overlap redshifts of $`15`$. However assuming feedback regulation of star formation in low mass galaxies, these histories require the unphysical value of $`f_{}f_{\mathrm{esc}}=1.0`$ in galaxies with $`v_{\mathrm{cir}}>v_{}`$. We therefore conclude that reionization by an early population of massive metal-free stars (with very large111We conservatively adopt the ionizing photon yield of $`100M_{}`$ metal-free stars. The yield is larger by a factor of $`2`$ for Pop-III stars with masses $`300M_{}`$ (Bromm, Kudritzki, & Loeb 2001). The use of the higher yield would have lead to the same $`\tau _{\mathrm{es}}`$ for values of $`f_{}f_{\mathrm{esc}}`$ at $`z>z_{\mathrm{tran}}`$ that are smaller by a factor of 2 than we quote. $`N_\gamma `$) must have contributed significantly to $`\tau _{\mathrm{es}}`$ at $`z15`$.

The formation of metal free stars is moderated with cosmic time by the gradual metal enrichment of the IGM. It has been argued that the dominant mode of star formation change to Pop-II once the metallicity of the associated gas rises above $`10^4`$Z (Schneider et al. 2002; Bromm et al. 2001). If all Pop-III stars formed in the mass range of $`140`$M$`{}_{}{}^{}260`$M where pair-unstable supernovae result in expulsion of metals into the surrounding gas (Heger & Woosley 2001), and if all metals were mixed uniformly throughout the entire IGM, then Pop-III star formation would have ceased as soon as $`0.3`$ ionizing photons were produced per cosmic hydrogen atom (Ciardi, Ferrara & White 2003). However, the highly inhomogeneous nature of the enrichment process (Furlanetto & Loeb 2003) and the formation of Pop-III stars beyond the above mass range, naturally lead to reionizion.

Must the Pop-III stars possess a heavy initial mass function (IMF)? Tumlinson & Shull (2000) computed the ionizing photon output per baryon of zero-metallicity stars with a Salpeter IMF and found only a 50% increase over the values for metal enriched Pop-II stars with the same IMF. This moderate increase is to be compared with the increase by more than an order of magnitude for the massive Pop-III stars described by Bromm, Kudritzki, & Loeb (2001) (the number of ionizing photons produced per baryon is nearly independent of stellar mass). The arguments presented above therefore imply that Pop-III stars were very massive rather than of a Salpeter IMF. This provides the first observational support to the results from hydrodynamic simulations of early star formation (Bromm et al. 2002; Abel et al. 2002).

In generating the star formation history we have extrapolated results from the local galaxy sample, assuming that feedback governs the observed stellar mass as a function of halo potential depth. However high redshift galaxies may have cooling times that are shorter than the dynamical time of their cold gas. This may lead to a starburst that might render the time-delayed supernovae feedback less effective. We now briefly consider this case. Assuming the star formation efficiency to be independent of halo mass (case A in Paper I), we find that reasonable values of $`f_{}f_{\mathrm{esc}}0.1`$ can yield $`\tau _{\mathrm{es}}0.17`$ in the absence of massive Pop-III stars. However, the corresponding reionization histories also have an early HII overlap. Therefore, unless the global product of escape fraction and star formation efficiency was more than an order of magnitude larger at $`z15`$ than at $`z9`$, our conclusion that an early generation of massive metal free stars is required to produce both $`\tau _{\mathrm{es}}0.17`$ and HII overlap at $`z9`$ is independent of whether feedback regulates star formation in low mass galaxies or not.

### 3.3 Reionization Histories

We now show that the joint constraints of $`\tau _{\mathrm{es}}=0.17\pm 0.04`$ and HII overlap at $`z9`$ can be met simultaneously for values of $`z_{\mathrm{tran}}>9`$ and $`f_{}f_{\mathrm{esc}}<1`$ if massive metal-free stars were present at high redshifts. In Figure 1 we show two cases of the evolution of $`Q`$ for hydrogen (solid lines), singly ionized helium (dashed lines) and doubly ionized helium (dotted lines). We assume222Salvaterra & Ferrara (2003) have recently presented the observational case that $`z_{\mathrm{tran}}9`$ through fitting the near-infrared background. The lower value of $`z_{\mathrm{tran}}`$ would allow a second HII overlap near $`z7`$ to be compatible with the large value of $`\tau _{\mathrm{es}}`$ and IGM temperature constraints. $`z_{\mathrm{tran}}=12`$ and $`f_{}f_{\mathrm{esc}}=0.08`$ at $`z<z_{\mathrm{tran}}`$, yielding the most recent HII overlap at $`z9`$. At $`z>z_{\mathrm{tran}}`$ we assume $`f_{}f_{\mathrm{esc}}=0.08`$ (left) and $`f_{}f_{\mathrm{esc}}=0.33`$ (right). The latter allows for the possibility of a higher star formation efficiency333Note that the star formation efficiency in the low-mass galaxies of interest, $`\eta _{}=0.02f_{}\left(M/10^8M_{}\right)^{2/3}[(1+z)/15]`$, is generally $`f_{}`$.  or escape fraction (Steidel et al. 2001) at high redshift. In the former case there is a significant reionization but failed HII overlap at early times, followed by a recombination and a successful HII overlap at $`z9`$ due to Pop-II stars. The two examples yield $`\tau _{\mathrm{es}}`$ within one standard deviation and at the most likely value of WMAP, respectively. Note that our model leads to the double reionization of helium by quasars at $`z3.5`$ independent of the contribution from Pop-III stars. This is in excellent agreement with current observational data \[e.g. Theuns et al. (2002) and references therein\]. The reionization of HeIII by quasars is characterized by a sharp rise in our model, with $`Q_{\mathrm{He}^{++}}`$ increasing from $`50\%`$ at $`z4.3`$ to $`100\%`$ at $`z3.7`$. We do not expect pre-heating from the reionization of HeIII to impact on the conclusion that the most recent HII overlap took place at $`z9`$.

Figure 1 also shows the evolution of $`Q`$ for hydrogen when helium is ignored (light solid line). In these calculations all photons blueward of 13.6eV reionize hydrogen. We find that the additional ionizing photons result in overlap of hydrogen at an earlier epoch, but the net optical depth to electron scattering is decreased.

So far we have assumed that atomic cooling dictates the mass of the smallest galaxies. However it has been suggested that significant production of X-rays could catalyze molecular hydrogen production allowing star formation in mini-halos before reionization (Haiman, Rees, & Loeb 1996; Cen 2002). The light dashed lines in the left hand panels of Figure 1 show the evolution in $`Q`$ for hydrogen (in the absence of helium) assuming $`T_{\mathrm{min}}=10^{2.4}`$K in the neutral IGM. We find that the feedback-regulation built into the model sufficiently limits the star formation efficiency in mini-halos so that they do not contribute significantly to the reionization history.

## 4 Conclusions

We have shown that Pop-II stars alone cannot reionize the universe sufficiently early to produce the value of $`\tau _{\mathrm{es}}17\pm 4\%`$ inferred from the WMAP data (Kogut et al. 2003; Spergel et al. 2003). This conclusion holds if feedback regulates star formation in early low-mass galaxies as observed in present-day dwarf galaxies. Without metal free stars, the condition that the most recent overlap of intergalactic HII regions must have occurred at $`z9`$ (Theuns et al. 2002) generically yields $`\tau _{\mathrm{es}}9\%`$. The joint requirements that $`\tau _{\mathrm{es}}17\%`$ and that HII regions overlapped at $`z9`$ imply that metal-free stars reionized a significant fraction of the IGM at $`z15`$. Moreover, we find that these stars must have been massive to provide the required order-of-magnitude increase in the production rate of ionizing photons relative to low mass Pop-II stars (Bromm, Kudritzki & Loeb 2001). Pop-III stars with a Salpeter IMF provide only a $`50\%`$ increase in $`N_\gamma `$ (Tumlinson & Shull 2000) and cannot account for the large $`\tau _{\mathrm{es}}`$.

For plausible reionization histories with metal-free stars we have found that singly and doubly ionized helium contributes as much as $`\mathrm{\Delta }\tau _{\mathrm{es}}=0.02`$, and should be included in reionization calculations. We also found that if the star formation efficiency is regulated in dwarf galaxies as observed in the local universe, then there is no significant contribution to the cosmic ionizing radiation field from halos with virial temperatures $`T_{\mathrm{vir}}10^4`$K, even if molecular hydrogen permits cooling of gas inside these mini-halos.

The evolution of the ionization state and metal enrichment of the IGM during the reionization epoch can be best probed by infrared spectroscopy of gamma-ray burst (GRB) afterglows (Furlanetto & Loeb 2002). These transient sources should be both abundant (Bromm & Loeb 2002) and detectable (Lamb & Reichart 2000; Ciardi & Loeb 2000) out to high redshifts. The Swift satellite, planned for launch this year, is expected to detect $`100`$ GRBs per year, out of which $`20\%`$ may originate at $`z5`$ (Bromm & Loeb 2002). Photometric identification of the Ly$`\alpha `$ trough can be used to select the high-redshift afterglows for which spectroscopic follow-up observations may be conducted on a 10-m class telescope (Loeb 2001).

After the submission of this paper, a number of different groups have also commented on the implications of the large $`\tau _{\mathrm{es}}`$ measured by WMAP. Some studies find that Pop-II stars are able to reionize the universe sufficiently early to produce $`\tau _{\mathrm{es}}0.13`$ (Ciardi, Ferrara & White 2003; Somerville & Livio 2003); however, this conclusion is challenged by the existence of the observed Gunn-Peterson trough in quasars at $`z6.3`$ (as noted by Ciardi et al. 2003) and by the temperature evolution of the IGM (Theuns et al. 2002; Hui & Haiman 2003). Other authors agree that the combined requirements of a large $`\tau _{\mathrm{es}}`$ and a late HII overlap implies that the efficiency of production of ionizing photons was higher at redshifts $`z10`$ (Haiman & Holder 2003; Sokasian, Abel, Hernquist & Springel 2003; Cen 2003). The conclusion that Pop-III stars needed to be massive to produce the observed $`\tau _{\mathrm{es}}`$ was also reached by Cen (2003) in agreement with our results.

We thank Rennan Barkana, Volker Bromm, Steve Furlanetto, Rachel Webster and the referee Andrea Ferrara for very useful comments on the manuscript. AL acknowledges generous support from the Institute for Advanced Study at Princeton and the John Simon Guggenheim Memorial Fellowship. This work was also supported in part by NSF grants AST-0071019, AST-0204514 and NASA grant ATP02-0004-0093.