text:
# 1 Introduction

## 1 Introduction

Phase transitions and coexistence of different phases in multi-component systems are of great importance to many areas of physics, chemistry and biology. An ideal system to study these phenomena is a multi-component dilute atomic gas Bose-Einstein condensate (BEC) mixture at zero temperature, due to the simplicity of its theoretical description. The mean-field approximation provides an excellent description of these systems. Other multi-component systems cannot be understood as well as these BEC mixtures, because their density is generally much higher, and the delta-function pseudopotential, used to describe interactions between atoms in BEC, is not appropriate for them. Instead, a true microscopic interaction potential must be employed to adequately describe such systems, hence modeling them is much harder.

Multi-component BECs have been extensively studied over the last few years -. These studies have been motivated by experimental work performed by the JILA  and MIT  groups. Many interesting effects have been experimentally determined and theoretically predicted, including topological properties of the ground and excited states, phase transitions and symmetry breaking , effects produced by a phase difference between components , stability properties , Josephson-type oscillations , four-wave mixing , and trapping of boson-fermion and fermion-fermion systems . Nevertheless, many features of BEC mixtures remain to be explored by theorists and experimentalists.

A large variety of different species can be used to produce mixtures of condensed bosons. Mixtures of two different elements, or of different isotopes of the same element, or simply different hyperfine states of the same atom  can be considered. Simulating experimental results for BEC binary mixtures requires knowledge of the scattering lengths of the atoms involved. To the extent that the values of the scattering lengths are known with insufficient accuracy, a full classification of different states is necessary within the range of possible values. This is also necessary in the context of tuning the scattering length, as can be done by means changing the external magnetic field near Feshbach resonances Ref. . Varying the scattering length, one can study phase transitions to states that break the symmetry of the trapping potential. Such states are known in the literature, and they were observed in the JILA experiment .

The classification of two-component condensates can also be carried out in systems with interconversion of components (“chemical reactions” between them), as in the case of an atom-molecule condensate, where the conserved quantity is the number of atoms plus twice the number of molecules (the number of atoms and number of molecules are not separately conserved). A mathematical model of the latter system can be formulated in terms of two coupled Gross-Pitaevskii equations (GPEs), which contain, in addition to the familiar cubic self- and cross-interaction nonlinear terms, quadratic terms that account for the “chemistry” (i.e., the interconversion) . In particular, an interesting prediction of the model is that a “soliton” state, i.e., a stationary self-supported condensate cloud (similar to “light bullets” in nonlinear optics ), may exist without any trapping potential present .

Many aspects of binary condensate mixtures have been treated in the literature, using various (mostly numerical) methods in order to predict results for various experimental setups (see, e.g., Refs. ,  and ). Nevertheless, a general classification of all the ground-state solutions is not yet available. This is understandable in view of many control parameters present in the models (three scattering lengths, particle numbers for both components and characteristics of the trap). A very general and elegant, but not explicit, algorithm for determining ground-state shapes has been proposed by Ho and Shenoy . Here we start with the same goal in mind, but also with the intention to provide a maximally straightforward and analytic set of expressions for the ground-state wavefunctions and energies of atomic gas BEC mixtures.

We use a computational method simulating the evolution of two-component GPEs in imaginary-time  in order to study phase separation of components in BEC mixtures. Results produced by this method are compared with analytical predictions based upon the Thomas-Fermi (TF) approximation applied to two-component GPEs. We analyze the changes in the structure of separated phase BEC mixtures with the variation of the $`s`$-wave scattering lengths and atom numbers. The changes can be predicted and understood using a simple TF approximation, which includes equations obtained from normalization conditions for both components. The TF picture is compared with solutions obtained using numerically simulated GPE evolution in imaginary time, which includes the kinetic energy of atoms neglected in the TF approximation. We find that, for the simple case of a spherically symmetric harmonic trapping potential, many spherically-symmetric phase-separated geometries are possible, depending on the ratios of the self- and cross- $`s`$-wave scattering lengths for atomic collisions. In the TF approximation, symmetry-broken phase-separated geometries (i.e., those whose symmetry is lower than that of the trap) are always energetically higher in energy than those with unbroken symmetry. Nevertheless, numerical simulations in imaginary-time show that a lower-symmetry state may be the lowest energy eigenstate, and thus determine a ground state of the system. Our method may be generalized to include a finite-range interaction between atoms. In the last section of the paper we study, by means of direct numerical simulations, how such interactions affect the geometry and shape of the ground state.

## 2 Mean-field description of two-component Bose-Einstein mixtures

In the present work, we concentrate on stationary states of BEC mixtures, (not their dynamics), therefore we start with the time-independent coupled GPEs, written in the standard notation:

$$\left(\mu _1\frac{\mathrm{}^2^2}{2m_1}+V_1(𝐫)+U_{11}|\psi _1(𝐫,t)|^2+U_{12}|\psi _2(𝐫,t)|^2\right)\psi _1(𝐫,t)=0,$$
(1)
$$\left(\mu _2\frac{\mathrm{}^2^2}{2m_2}+V_2(𝐫)+U_{12}|\psi _1(𝐫,t)|^2+U_{22}|\psi _2(𝐫,t)|^2\right)\psi _2(𝐫,t)=0.$$
(2)
Here $`\mu _{1,2}`$ are chemical potentials of the two species, and $`V_{1,2}(𝐫)`$ are two isotropic parabolic trapping potentials, i.e.,

$$V_j(𝐫)=(m_j/2)\omega _j^2r^2,j=1,2,$$
(3)
$`\omega _j`$ being the corresponding frequencies of harmonic oscillations of a trapped particle . Further, $`U_{ij}\left(4\pi \mathrm{}^2/m_{ij}\right)a_{ij}`$ are atom-atom interaction strengths, proportional to the $`s`$-wave scattering lengths $`a_{11}`$, $`a_{22}`$, and $`a_{12}`$ for the $`1+1`$, $`2+2`$, and $`1+2`$ collisions, respectively, where $`1`$ and $`2`$ numerate the components, and $`m_{ij}=m_i`$ if $`i=j`$ and $`m_{ij}=m_1m_2/\left(m_1+m_2\right)`$ if $`ij`$. For simplicity, in the numerical calculations presented here we take $`m_1=m_2m`$, and assume that the magnetic moments of atoms belonging to the different components are equal, so that the corresponding trapping potentials are equal too, $`\omega _1=\omega _2\omega `$, but this condition as well as the spherical symmetry condition may be readily relaxed by means of rescaling variables. Furthermore, we assume that the scattering lengths are real, i.e., we assume that collisions are not lossy. We also assume that all the scattering lengths for both different and alike atoms are positive; otherwise, the classification of the possible states becomes very cumbersome.

Our calculations were carried out, simulating the evolution in the time-dependent GPEs in imaginary time , so that to let the solution relax to the ground state. The computations used the split operator method with the fast Fourier transform, similar to that used in Ref. . The chemical potentials are obtained by computing the net energy (the sum of kinetic, potential, self- and cross- nonlinear mean-field energies) for each component. We have chosen the wavefunctions $`\psi _1(𝐫,t)`$ and $`\psi _2(𝐫,t)`$ to be normalized to the number of particles in each component, so that $`|\psi _i(𝐫,t)|^2d^3𝐫=N_i`$. Choosing the symmetry of an initial configuration in the imaginary-time simulations, a solution $`\psi _{1,2}(𝐫)`$ which minimizes the total energy

$$=d^3𝐫[\mu _1|\psi _1(𝐫,t)|^2+\mu _2|\psi _2(𝐫,t)|^2\frac{U_{11}}{2}|\psi _1(𝐫,t)|^4\frac{U_{22}}{2}|\psi _2(𝐫,t)|^4U_{12}|\psi _1(𝐫,t)|^2|\psi _2(𝐫,t)|^2],$$
(4)
within the class of functions possessing this symmetry can be found. The ground state of the two-component Hamiltonian is the one with the smallest value of $``$ for different symmetry classes.

## 3 Thomas-Fermi Approximation

The TF approximation can be used to describe a zero temperature condensate in the cases when the trapping-potential and mean-field nonlinear terms in GPEs are attractive and repulsive respectively, and the number of atoms is large so that the mean-field energies are large compared to the kinetic energy. For the two-component system with overlapping wavefunctions,

$$\left(\begin{array}{c}\mu _1V_1(𝐫)\\ \mu _2V_2(𝐫)\end{array}\right)=\left(\begin{array}{cc}U_{11}& U_{12}\\ U_{21}& U_{22}\end{array}\right)\left(\begin{array}{c}|\psi _1(𝐫,t)|^2\\ |\psi _2(𝐫,t)|^2\end{array}\right).$$
(5)
These equations can be solved as a linear system of equations for $`|\psi _1(𝐫,t)|^2`$ and $`|\psi _2(𝐫,t)|^2`$ in terms of the chemical potentials $`\mu _1`$ and $`\mu _2`$ to obtain:

$$|\psi _1(𝐫)|^2=\frac{\left[\mu _1U_{11}\mu _2U_{12}\right]\left[U_{22}U_{12}\right](m/2)\omega ^2r^2}{U_{11}U_{22}U_{12}U_{12}},$$
(6)
$$|\psi _2(𝐫)|^2=\frac{\left[\mu _2U_{11}\mu _1U_{12}\right]\left[U_{11}U_{12}\right](m/2)\omega ^2r^2}{U_{11}U_{22}U_{12}U_{12}}.$$
(7)
The chemical potentials $`\mu _1`$ and $`\mu _2`$ are determined from the normalization conditions, $`N_i=d^Dx|\psi _i(𝐫,t)|^2`$, where $`D`$ is the dimension. We shall plot examples for $`D=1`$, but our numerical method is valid for higher dimensions as well, hence we derive all the formulas for the general case.

When phase separation occurs and there are regions in the physical space occupied by one component only, the TF approximation leads, instead of Eqs. (5), to the corresponding one-component TF equations. For example, if a phase with only species $`1`$ exists in a particular region of space, an equation of the form

$$|\psi _1(𝐫,t)|^2=(\mu _1V_1(𝐫))/U_{11},$$
(8)
is to be used in this region. If another phase exists wherein the species $`1`$ and $`2`$ are mixed, Eqs. (5) are relevant for that region. The chemical potentials $`\mu _1`$ and $`\mu _2`$ must be determined by setting the number of atoms of each type equal to the integral of the corresponding density over the whole space.

Inspection of Eqs. (5) suggests then that all the solutions can be classified according the signs of three parameters: $`detUU_{11}U_{22}\left(U_{12}\right)^2`$, and $`\alpha _jU_{jj}U_{12}`$. In particular, $`\alpha _{1,2}`$ determine signs of the curvature (coefficients in front of $`r^2`$) of the effective quadratic potentials for the two components in Eqs. (6) and (7). Qualitatively different types of possible states with overlapping wave functions (i.e., disregarding regions where only one of the species is present) identified by the TF analysis are defined in Table 1.

Let us focus on those cases when kinetic energy contribution does not change the general structure of the solution, but only generates narrow transient layers on the scale of the corresponding healing length, as in the single-component case when the TF approximation is valid (thus we consider large-size condensates. The analysis will include the TF configurations with and without the overlap of the two different condensate wavefunctions. As already mentioned, in the most cases both wavefunctions do not overlap everywhere, i.e., there is a region where only one wavefunction is different from zero. Consequently, search for the lowest-energy state of the mixture cannot rely solely on Eqs. (6)-(7) obtained in the assumption that the overlapping takes place everywhere.

Combining the cases represented in Table 1 and single-wavefunction solutions within the TF approximation, we distinguish two general types of solutions: unseparated ones, having an overlap region where both wavefunctions are nonzero, and separated solutions which do not contain any overlap region. In the latter case, we shall see from the analysis of a full GP equation (with kinetic energy included) and also a nonlocal version of the two-species model, with a finite range of the interatomic interactions, that it is necessary to further distinguish between weak ($`U_{11}U_{22}\left(U_{12}\right)^2`$) and strong ($`U_{11}U_{22}\left(U_{12}\right)^2`$) separation .

### 3.1 Partially Overlapping Wavefunctions $`(U_{11}U_{22}\left(U_{12}\right)^2>0)`$

In the case $`detU>0`$ (cases 1, 3 and 4 in Table 1), we have checked numerically that the minimum-energy solution is given by the wavefunctions of the form shown in Fig. 1. Near the origin, both wavefunctions coexist up to the point where one of them vanishes. Past this point, one wavefunction vanishes, while the other one remains nonzero, following the single-component solution, Eq. (8). Fig. 1 shows two different cases that are possible with the scenario described above. In Fig. 1a, the two effective trapping potentials have the same sign of their curvature, i.e., $`\alpha _1\alpha _2>0`$, in the overlap region (this is case 1 in Table 1). Fig. 1b presents another situation, when the two effective potentials have opposite curvatures, $`\alpha _1\alpha _2<0`$ (these are cases 3 and 4 in Table 1).

### 3.2 Separated Wavefunctions $`(U_{11}U_{22}\left(U_{12}\right)^2<0)`$

This category is represented by cases 2, 5 and 6 from Table 1. The simplest configuration is that with one wavefunction being different from zero in the region around the origin and vanishing beyond a separation radius, $`R`$, while the second component surrounds the first one. In this case, we can express $`\mu _1`$ and $`\mu _2`$ as functions of $`R`$ and minimize the net free energy $``$, in order to find the lowest eigenstate of this type. The normalization conditions for the wavefunctions of the two components give a set of relations between $`R`$, the chemical potentials $`\mu _1`$ and $`\mu _2`$, and the number of atoms in each condensate:

$$N_1=_0^Rd^Dr\left[\mu _1V(𝐫)\right]/U_{11},$$
(9)
$$N_2=_R^{R_0}d^Dr\left[\mu _2V(𝐫)\right]/U_{22}.$$
(10)
Here $`V(𝐫)`$ is the binding potential (3), and $`R_0`$ is a outer radius at which the wavefunction of the second component vanishes in the FT approximation. We first consider the 1D case and then show how these considerations can be generalized to two and three dimensions.

#### 3.2.1 One-Dimensional Case

To find the value of the radius $`R`$ minimizing $``$, we solve the set of the coupled equations (9) and (10) for $`\mu _1`$ and $`\mu _2`$. The first equation can be solved directly to yield $`\mu _1`$ as a function of $`R`$. The second equation is more complicated – it can be solved analytically only in the 1D and 2D cases, but not in 3D. In the 1D (3D) case, one needs to solve a third- (fifth-) order algebraic equation to find $`\mu _2`$ as a function of $`R`$. The system of equations (9) and (10) in 1D reduces to

$$2(\mu _1R\frac{1}{6}m\omega ^2R^3)=U_{11}N_1$$
(11)
$$2(\frac{(2\mu _2)^{3/2}}{3\sqrt{m\omega ^2}}\mu _2R+\frac{1}{6}m\omega ^2R^3)=U_{22}N_2.$$
(12)
In the case under consideration here, without spatial overlap of components, the total energy is simply a sum of the average values of harmonic potential and half of the nonlinear term in the GP equation for each component state: $`=\frac{1}{2}_i<\psi _i|m\omega ^2x^2+U_{ii}|\psi _i|^2/2|\psi _i>`$. If we substitute the direct expression for the wavefunction in the TF approximation we obtain:

$$=(\frac{(2\mu _2)^{5/2}}{5\sqrt{m\omega ^2}}\mu _2^2R+\frac{1}{20}(m\omega ^2)^2R^5)/U_{22}+(\mu _1^2R+\frac{1}{20}(m\omega ^2)^2R^5)/U_{11}$$
(13)
and the chemical potentials can be found using Eqs. (11)-(12).

#### 3.2.2 Generalization to Two and Three Dimensions

A particularly simple form of these equations is obtained upon introducing a TF radius $`\left(R_{\mathrm{TF}}\right)_{1,2}`$ of the condensate in the corresponding dimension. The TF radius is defined as a radius of the single spherically-symmetric condensate obtained in the TF approximation. We can find an explicit dependence between the chemical potentials of the condensates and the separation radius of the two phases in 2D. In this case, we again define the TF radii, which in 2D are equal to: $`\left(R_{\mathrm{TFi}}\right)^4=8N_iU_{ii}/(\pi m\omega ^2)`$. In this case, we obtain the following set of equations for chemical potentials and total energy:

$$\mu _1=\frac{m\omega ^2}{4R^2}\left[\left(R_{\mathrm{TF1}}\right)^4+R^4\right],$$
(14)
$$\mu _2=\frac{m\omega ^2}{2}\left[\left(R_{\mathrm{TF2}}\right)^2+R^2\right],$$
(15)
$$=\frac{1}{(m\omega ^2)^2}(\frac{2\mu _1^2r^2\frac{1}{6}r^6}{R_{TF1}^4}+\frac{\frac{8}{3}\mu _2^32\mu _2^2r^2+\frac{1}{6}r^6}{R_{TF2}^4}).$$
(16)
where r is defined as $`r=R\sqrt{m\omega ^2}`$.

In 3D we can define the TF radius is given by $`\left(R_{\mathrm{TFi}}\right)5=(15U_{ii})/(2\pi m\omega ^2)`$, and we can derive a set of equations which can be solved for the chemical potentials vs. the separation radius:

$$\mu _1=\frac{m\omega ^2}{10R^3}\left[3R^5+2\left(R_{\mathrm{TF1}}\right)^5\right],$$
(17)
$$4(\frac{2\mu _2}{m\omega ^2})^{5/2}+3R^55\frac{2\mu _2}{m\omega ^2}R^32\left(R_{\mathrm{TF2}}\right)^5=0.$$
(18)
The energy is given in terms of the separation radius R by:

$$=\frac{1}{(m\omega ^2)^{5/2}}\left(\frac{\frac{5}{2}\mu _1^2r^3\frac{15}{56}r^7}{R_{TF1}^5}+\frac{\frac{5(2\mu _2)^{7/2}}{14}\frac{5}{2}\mu _2^2r^3+\frac{15}{56}r^7}{R_{TF2}^5}\right)$$
(19)
where $`r=R\sqrt{m\omega ^2}`$.

#### 3.2.3 Symmetry Breaking Solutions

Eigenstates of the binary-condensate system that break the symmetry of the trapping potential exist. In 1D, a solution of this kind is given by TF parabolas that are stuck together. An example is shown in Fig. 2. In this case one can derive equations for $`\mu _1`$ and $`\mu _2`$ in the same way as in Sec. 3.2.1, integrating the densities and substituting the result into the normalization conditions. The generalization to the dimensions higher than one may be only obtained if the separation surface (which reduces in 1D to a single point) is simple. In 1D, the corresponding coupled equations take the form

$$\frac{(2\mu _1)^{3/2}}{3\sqrt{m\omega ^2}}\mu _1R+\frac{1}{6}m\omega ^2R^3=U_{11}N_1$$
(20)
$$\frac{(2\mu _2)^{3/2}}{3\sqrt{m\omega ^2}}+\mu _2R\frac{1}{6}m\omega ^2R^3=U_{22}N_2$$
(21)
$$=(\frac{(2\mu _1)^{5/2}}{5\sqrt{m\omega ^2}}\mu _1^2R+\frac{1}{20}(m\omega ^2)^2R^5)/U_{11}+(\frac{(2\mu _1)^{5/2}}{5\sqrt{m\omega ^2}}+\mu _2^2R\frac{1}{20}(m\omega ^2)^2R^5)/U_{22}$$
(22)
We have checked numerically that this solution cannot give rise to a minimum of the free energy, hence, within the framework of the FT approximation, the ground state cannot be the one with broken symmetry. The difference in the energy between symmetric and asymmetric cases is usually very small making them almost degenerate. The degeneracy is exact in the limit when $`U_{11}U_{22}`$ and is removed by the kinetic energy.

## 4 The role of kinetic energy

In this section we present results of our studies of the contribution of the kinetic energy on the total energy and on the functional dependence of the ground state of the binary mixtures of the BEC. For a single condensate, in the regime of the validity of TF approximation, the kinetic energy creates a healing length, $`\xi `$, in the region where the condensate wavefunction tends to zero. This healing length is of the order of $`(8\pi na)^{1/2}`$ , where $`a`$ is a scattering length and $`n`$ is an average density of the condensate. For the binary mixtures of BECs there are two length scales that can be defined in order to characterize two kinds of boundary regions. One is the ordinary healing length of single condensate and refers to the healing of the wavefunction outside region of the coupled wavefunctions. But for mixtures another region between the two components exist and a penetration depth, $`\chi `$, as a length scale over which two components overlap. The penetration depth is a function of $`detU`$. For $`detU<0`$ the lowest energy state consists of partially overlapping wavefunctions; hence the penetration depth is of order of the size of the condensate. With decreasing $`detU`$ the penetration depth becomes smaller and goes to zero in the limit of strong repulsion as $`detU\mathrm{}`$. At the same time, the contribution of the kinetic energy to the total energy becomes more important, in spite of the shrinking overlap region. This situation is illustrated in Fig. 3. The energy of the lowest eigenstate is plotted in the symmetric and asymmetric classes as a function of $`detU/(U_{11}U_{22})`$. Fig. 3 Only one curve is plotted for $`detU>0`$ where the contribution of the kinetic energy is negligible and Thomas-Fermi approximation gives an excellent prediction for both the value of the ground state energy and its wavefunction. As the value of $`detU`$ becomes negative, the lowest energy state within a TF approximation consists of two separated components and the energy does not depend on $`U_{12}`$ and is plotted as a horizontal dashed line in Fig. 3. $`detU<0`$ the contribution of the kinetic energy is substantial and it is larger for the symmetric case which has two interfaces between phases. The asymmetric configuration has only one interface in the lowest energy state. Hence, the ground state looses symmetry of the trap.

In order to search for stable symmetry-breaking solutions we kept $`N_1`$,$`U_{11}`$ and $`U_{12}`$ constant and varied $`N_2`$ and $`U_{22}`$. Fig. 4 plots the ratio of the total energy in the asymmetric case to the energy of the symmetric one as a function of these two variables. Almost all solutions are symmetry-breaking ones and a trough is formed near $`U_{22}=U_{11}`$. The trough is an optimal region for finding symmetry-breaking solutions.

## 5 Finite Interaction Range

We have so far considered the mean-field description of a BEC mixture assuming a zero-range delta-function pseudopotential. It is of interest to consider the effects of a finite-range interatomic interaction on the ground-state structure in the two-component condensate. We introduce a pseudopotential in the form of a normalized Gaussian with a finite width (i.e., range) which recovers the zero-range limit result as the range vanishes. We search for changes in the structure of the ground state of the two-component system as the range of the intercomponent interactions only are varied, keeping the delta-function pseudopotential for the self-interactions. The results displayed below were obtained by means of direct numerical simulations, not the FT approximation. In 1D, the intercomponent-interaction terms in Eqs. (1) $`U_{12}|\psi _i(x,t)|^2`$ are replaced by a nonlocal expression,

$$U_{12}\frac{1}{\sqrt{2\pi d^2}}_{\mathrm{}}^+\mathrm{}𝑑y\mathrm{exp}(\frac{(xy)^2}{2d^2})|\psi _i(y,t)|^2,$$
(23)
where $`d`$ is the interaction range. The Gaussian form was chosen to model a finite-range potential solely for its simplicity (see also Ref.  for the role of a finite interaction range in attractive single-component BEC within this model). This family of finite-range potentials has a constant scattering length within the first Born approximation. All the cases that we investigate below correspond to configurations with separated wavefunctions in the usual TF limit (see Sec. 3.2).

### 5.1 Weak separation $`(U_{11}U_{22}\left(U_{12}\right)^2)`$

Here we discuss case 2 of Table 1. Fig. 5 shows the overlap region between the wavefunctions of the two components; this overlap region grows as the interaction range increases. Starting from relatively well separated phases, we end up with a complete overlap of the two components (i.e., the component located initially outside the narrow-width component finally penetrates it). The intercomponent interaction parameters are chosen as $`U_{12}=1.02U_{11}=1.02U_{22}`$, which places this case not far from the boundary between the cases of separated and overlapping phases. In other words, increasing the interaction range forces a transition from separated phases to penetrating ones. Within the zero-range model, this would correspond to a transition from $`U_{12}\sqrt{U_{11}U_{22}}`$ to $`U_{12}\sqrt{U_{11}U_{22}}`$, i.e., to attenuation of the interaction. A simple argument justifies this conclusion. Suppose the interaction range is so large that the long-range potential varies only slightly over the extent of the interface between the nearly separated components (where the intercomponent interaction is important). Then, the interaction terms of Eq. (23) would be approximately constant, producing only a shift in the total energy of the mixture, but not affecting the shape of the two wavefunctions. When $`d`$ becomes comparable to the penetration depth (see also Ref. ), phase separation is reduced. This is clearly illustrated in Fig. 5.

Although the parameters in the 1D solution are not directly relevant to an experimental situation, we decreased the range of the intercomponent interactions so that they are in a realistic range of values. From the aforementioned arguments, we see that in order to observe a difference from the zero-range case, the interaction range should be comparable to the penetration depth. In order for the interaction range to be small and yet correspond to the boundary between overlapping and segregated phases, a large disparity in the number of atoms in the two components is required. This situation is depicted in Fig. 6, where $`U_{12}=1.01U_{11}=1.01U_{22}`$, but $`N_1N_2`$ . We observe a qualitative change in the ground-state solution: the component that initially was at the center of the trap moves to its outskirts as the interaction range grows. Although $`d=0.1\sqrt{\mathrm{}/m\omega }`$ (bottom frame in Fig. 6) would usually correspond to several tens of nanometers, we argue that one can optimize the sensitivity to $`d`$ by considering the regime of parameters near the boundary between the penetrating and segregated phases. The possibility of manipulating the strength of atomic collisions is not excluded, as Feshbach resonances have been observed in BEC samples , and several other proposals in this respect have been put forward . Using such techniques one could vary the intercomponent scattering length in order to scan the region near $`detU=0`$. Then, by comparing the measured structure of the mixture to the predictions of the two theoretical models (i.e., the ones with zero and finite interaction range) one could determine the effective range of interactions (which is a parameter in the latter model). Thus, one can probe the microscopic parameter $`d`$ via a magnified effect such as the qualitative change of the condensate structure from to separated phase to penetrating phase.

### 5.2 Strong separation $`(U_{11}U_{22}\left(U_{12}\right)^2)`$

In the preceding section we demonstrated the effective attenuation of the mean-field repulsive interaction between two components of a BEC due to an increase in the range of the interaction for the case when $`U_{11}U_{22}(U_{12})^2`$. Now we turn to the case when the parameters of the BEC mixture are far from threshold for the onset of the penetrating phase. In Fig. 7 is for parameters $`U_{12}=7U_{11}=7U_{22}`$, hence the interface between the two components is very sharp. Since now it is much easier to match the range of interactions with the penetration depth, one might expect that effects similar to those described above will appear at even smaller values of $`d`$. However, this is not the case. The mutual repulsion of the components remains very strong even if reduced by a finite interaction range. As the interaction range increases, the two components tend to move apart, yielding two completely separated phases.

## 6 Summary and Conclusion

We have presented a detailed classification of stable solutions for binary mixtures of dilute atomic condensates. The analysis is particularly simple within the Thomas-Fermi approximation. Within this approximation one can distinguish two general classes of ground state for two component condensate mixtures: unseparated ones with an overlap region (both component wavefunctions are simultaneously nonzero) and separated ones not containing an overlap region (except for the tail penetration). The latter contains also solutions that do not posses the symmetry of the trapping potential. Components are separated if $`detUU_{11}U_{22}<U_{12}^2<0`$ and they overlap if $`detU>0`$. The predictions from the TF approximation become ambiguous in the region of parameters where phase separated solutions which break the symmetry of the trap are energy-degenerate with the phase-separated solutions preserving the symmetry. In this case, it is crucial to include the contribution of the kinetic energy operator and of the mutual interaction energy in determining the structure of the ground state geometry, which tend to favour the asymmetric solutions. The physical reason for this is probably a smaller interface region in the asymmetric solutions. It is these interface regions which contribute most to both the kinetic and mutual interaction energies. We have carried out our numerical calculations in 1D, but our conclusions should be valid in two and three dimensions as well. The condensates, if not overlapping, should have a propensity towards states with minimal interface surface area. In the $`detU<0`$ case, with kinetic energy included, we further recognize a weak separation regime ($`detU0`$) and strong separation regime ($`detU0`$) when a penetration depth goes to zero. Contribution of the kinetic energy to the total energy increases with decreasing $`detU`$ (for negative $`detU`$), in spite of the decreasing interface size (see Fig. (3). This is due to the increasing importance of the $`U_{12}`$ term that gives the cross-interaction energy of atoms from different components.

Since the size of the overlap region can be very small (smaller than a single condensate healing length), it is of interest to investigate the possible impact of a small, but non-zero, interaction range in a binary condensate. We have developed a model of finite range potential by introducing a pseudopotential in the form of a normalized Gaussian with finite width. We identified two distinctly different cases. In one case, that of significant overlap, the finite range tends to increase the penetration over the delta function interaction, and in the other, that of strong separation, the finite range leads to a trough between the two condensates.

We thank M. Gajda and J. Mostowski for stimulating discussions. This work was supported in part by the US-Israel Binational Science Foundation and the James Franck Binational German-Israel Program in Laser-Matter Interaction (Y.B.B.). K.R. and K.G. acknowledge the support of the subsidy from the Foundation for Polish Science and of the Polish KBN Grant no 2 P03B 057 15. Part of the results have been obtained using computers at the Interdisciplinary Center for Mathematical and Computational Modeling (ICM) at Warsaw University.


format: text

text:
# Prompt Photons in Photoproduction at HERA

## 1 Introduction

Isolated high transverse energy (“prompt”) photon processes at HERA could yield information about the quark and gluon content of the photon, together with the gluon structure of the proton. The particular virtue of prompt photon processes is that the observed final state photon emerges directly from a QCD diagram without the subsequent hadronisation which complicates the study of high $`E_T`$ quarks and gluons.

In a ZEUS paper  the observation of prompt photons was first confirmed at HERA. More recently , ZEUS collaboration has measured the cross sections of inclusive prompt photons in photoproduction reactions, using an integrated luminosity of 38.4 pb<sup>-1</sup>. Comparisons are made with predictions from Monte Carlo models containing leading-logarithm parton showers, and with next-to-leading-order QCD calculations, using currently available parameterisations of the photon structure.

## 2 Evaluation of the photon signal

The data used here were obtained from $`e^+p`$ running in 1996–97 at HERA, with $`E_e=27.5`$ GeV, $`E_p=820`$ GeV.

The major components in the analysis are the central tracking detector(CTD) and the uranium calorimeter(UCAL). Prompt photons are detected in the barrel section of the calorimeter, which consists of an electromagnetic section (BEMC) followed by two hadronic sections. It enable a partial discrimination between single $`\gamma `$ signals and the decay product of neutral mesons. A typical high-$`E_T`$ photon signal is observed in a small cluster of BEMC cells, with no associated CTD track. An isolation cone was also imposed around photon candidates within a cone of unit radius in $`(\eta ,\varphi )`$, to reduce backgrounds from dijet events with part of one jet misidentified as a single photon.

Two shape-dependent quantities were studied in order to further distinguish $`\gamma `$, $`\pi ^0`$ and $`\eta `$ signals. These were (1) the mean width $`<\delta Z>`$ of the BEMC cluster in $`Z`$ and (2) the fraction $`f_{max}`$ of the cluster energy found in the most energetic cell in the cluster. The $`<\delta Z>`$ distribution is shown in Figure 1(a), in which peaks due to the $`\gamma `$ and $`\pi ^0`$ contributions are clearly visible. The tail quantified the $`\eta `$ background; photon candidates in this region were removed.

The extraction of the photon signal from the mixture of photons and a neutral meson background was done by means of the $`f_{max}`$ distribution. Figure 1(b) shows the shape of the $`f_{max}`$ distribution for the final event sample, after the $`<\delta Z>`$ cut, fitted to the $`\eta `$ component determined from the $`<\delta Z>`$ distribution and freely-varying $`\gamma `$ and $`\pi ^0`$ contributions. Above an $`f_{max}`$ value of 0.75, the distribution is dominated by the photons; below this value it consists mainly of meson background. The numbers of candidates with $`f_{max}0.75`$ and $`f_{max}<0.75`$ were calculated for the sample of events occurring in each bin of any measured quantity. From these numbers, and the ratios of the corresponding numbers for the $`f_{max}`$ distributions of the single particle samples, the number of photon events in the given bin was evaluated. Further details of the background subtraction method are given in reference.

## 3 Results

We evaluate cross sections for prompt photon production corrected by means of PYTHIA using GRV photon structure functions . A bin-by-bin correction factors were applied to the detector–level measurements so as to correct to cross sections in the $`\gamma p`$ centre–of–mass energy 134 – 285 GeV. The systematic error of 15% were taken into account and were finally combined in quadrature. The main contributions are from the energy scale of the calorimeter and the background subtraction. In presenting cross sections, comparison is made with two types of theoretical calculation, in which the pdf sets taken for both the photon and proton can be varied. These are (1) PYTHIA and HERWIG calculations evaluated at the final-state hadron level and (2) NLO parton–level calculations of Gordon (LG) and of Krawczyk and Zembrzuski (K&Z).

Figure 2 (a) gives the inclusive cross–section $`d\sigma /dE_T^\gamma `$ for isolated prompt photons in the range $`0.7<\eta ^\gamma <0.9`$. All the theoretical models describe the shape of the data well; however the predictions of PYTHIA and especially HERWIG are too low in magnitude. The LG and K&Z calculations give better agreement with the data.

The inclusive cross–section $`d\sigma /d\eta ^\gamma `$ for isolated prompt photons in the range $`5<E_T^\gamma <10`$ GeV is shown in figure 2 (b) and compared with theoretical calculations, using two sets of photon pdf, GS and GRV. The LG and K&Z calculations gives a good description of the data for forward (proton direction) $`\eta ^\gamma `$ range and are similar to PYTHIA prediction. However all the calculations lie below the data in the lower $`\eta ^\gamma `$ range, where the curves using the GS parton densities give poorer agreement than those using GRV.

The discrepancy between data and theory at negative $`\eta ^\gamma `$ is found to be relatively strongest at low $`\gamma p`$ centre-of-mass energy range. In the lowest $`W`$ range (134–170 GeV), both theory and data show a peaking at negative $`\eta ^\gamma `$, but it is stronger in the data. In the highest $`W`$ range (212–285 GeV), agreement is found between theory and data. The movement of the peak can be qualitatively understood by noting that for fixed values of $`E_T`$ and $`x_\gamma `$, where $`x_\gamma `$ is the fraction of the incident photon energy that contributes to the resolved QCD subprocesses, measurements at increasing $`y`$ correspond on average to decreasing values of pseudorapidity. By varying the theoretical parameters, the discrepancy was found to correspond in the K&Z calculation to insufficient high $`x_\gamma `$ partons in the resolved photon.

## 4 Conclusions

The photoproduction of isolated prompt photons within the $`\gamma p`$ centre-of-mass energy range 134–285 GeV has been measured in the ZEUS detector at HERA. Inclusive cross sections have been presented as a function of $`E_T^\gamma `$ for photons in $`0.7<\eta ^\gamma <0.9`$, and as a function of $`\eta ^\gamma `$ for photons with $`5<E_T^\gamma <10`$ GeV.

Comparisons have been made with predictions from LO Monte Carlos, and from NLO calculations. The models are able to describe the data well for forward $`\eta ^\gamma `$, but are low in the rear direction. None of the available variations of the model parameters was found to be capable of removing the discrepancy with the data. This result would appear to indicate a need to review the present theoretical modelling of the parton structure of the photon.

## Acknowledgments

We are grateful to L. E. Gordon, Maria Krawczyk and Andrzej Zembrzuski for helpful conversations, and for providing theoretical calculations.


text:
# 1 Higgs boson decay rates in GeV for different Higgs masses in GeV.

## Acknowledgements

We should like to thank G. Duckeck, B. Kniehl, A. Leike, A. Martin, M. Stanitzki and P. Zerwas for useful discussion.


text:
# 1 Introduction

## 1 Introduction

Measurements of the individual $`b`$ hadron lifetimes represent an important test of the present knowledge of nonspectator effects in the $`b`$ hadron decay dynamics, such as Pauli interference, $`W`$ exchange and weak annihilation. Based on the heavy quark expansion formalism, the difference between the lifetimes of the $`b`$ baryons and mesons is predicted to depend on $`1/m_b^2`$ and higher order terms, whereas meson-meson differences depend only on $`1/m_b^3`$ and higher order terms . The predicted hierarchy is $`\tau _{\mathrm{\Lambda }_b}<\tau _{\overline{B}^0}\tau _{B_s}<\tau _B^{}`$. Differences are expected to be at the level of a few percent, which sets the scale of the experimental precision required.

This paper reports an improved measurement of the $`\overline{B}^0`$ and $`B^{}`$ lifetimes with the Aleph detector at LEP, using approximately four million hadronic decays of the Z, collected in the period 1991–1995. This data sample was recently reprocessed, achieving higher efficiency and better resolution in the track reconstruction, which is highly beneficial for this analysis.

Semileptonic decays of $`\overline{B}^0`$ and $`B^{}`$ mesons are partially reconstructed by identifying events containing a lepton (electron or muon) with an associated $`D^0`$ or $`D^+`$ meson. The resulting $`D^0`$-lepton ($`D^0\mathrm{}^{}`$) and $`D^+`$-lepton ($`D^+\mathrm{}^{}`$) event samples consist mostly of $`B^{}`$ and $`\overline{B}^0`$ decays, respectively (charge conjugate modes are implied throughout this paper).

Previous measurements of the $`\overline{B}^0`$ and $`B^{}`$ lifetimes are reported in .

## 2 The ALEPH detector

A detailed description of the Aleph detector can be found elsewhere . A high resolution vertex detector (VDET) consisting of two layers of silicon with double-sided readout provides measurements in the $`r\varphi `$ and $`z`$ directions at average radii of 6.5 cm and 11.3 cm, with $`12`$ $`\mu `$m precision at normal incidence. The VDET provides full azimuthal coverage, and polar angle coverage to $`|\mathrm{cos}\theta |<0.85`$ for the inner layer only and $`|\mathrm{cos}\theta |<0.69`$ for both layers. Outside the VDET, particles traverse the inner tracking chamber (ITC) and the time projection chamber (TPC). The ITC is a cylindrical drift chamber with eight axial wire layers at radii of 16 to 26 cm. The TPC measures up to 21 space points per track at radii between 40 and 171 cm, and also provides a measurement of the specific ionization energy loss ($`dE/dx`$) of each charged track. These three detectors form the tracking system, which is immersed in a 1.5 T axial magnetic field provided by a superconducting solenoid. The combined tracking system yields a momentum resolution transverse to the beam axis of $`\sigma (p_T)/p_T=6\times 10^4p_T0.005`$ ($`p_T`$ in GeV/$`c`$). The resolution of the three-dimensional impact parameter for tracks having two VDET hits can be parametrized as $`\sigma =25\mu \text{m}+95\mu \text{m}/p`$ ($`p`$ in GeV/$`c`$.

The electromagnetic calorimeter (ECAL) is a lead/wire-chamber sandwich operated in proportional mode. The calorimeter is read out in projective towers that subtend typically $`0.9^{}\times 0.9^{}`$ in solid angle, segmented in three longitudinal sections. The hadron calorimeter (HCAL) uses the iron return yoke as absorber. Hadronic showers are sampled by 23 planes of streamer tubes, with analogue projective tower and digital hit pattern readout. The HCAL is used in combination with two double layers of muon chambers outside the magnet for muon identification.

Recently the LEP1 data were reprocessed using improved reconstruction algorithms. The features that are particularly relevant for the enhancement of the charmed meson reconstruction efficiency are the following. A new VDET pattern recognition algorithm allows groups of several nearby tracks which may share common hits to be analysed together, to find the hit assignments that minimize the overall $`\chi ^2`$ for the event. The improvement in the hit association efficiency is more than 2% (from 89.2% to 91.0% in $`r\varphi `$ and from 85.6% to 88.2% in $`z`$). Information on the drift time from the TPC wires is combined with that obtained from the pads to reduce the error in the $`z`$ coordinate by a factor of two. A 30% improvement in the $`r\varphi `$ coordinate resolution is achieved for low momentum tracks by correcting the pad coordinates for ionisation fluctuations along the tracks as measured by the wires. The particle identification ($`dE/dx`$) is improved by combining pulse height data from the TPC pads with that of the wires. The improvements in the charmed meson reconstruction efficiencies with respect to the previous analysis range from 10 to 30%, depending on the decay channel.

## 3 Event selection

The $`D^+\mathrm{}^{}`$ and $`D^0\mathrm{}^{}`$ event samples consist of an identified lepton ($`e`$ or $`\mu `$) associated with a fully reconstructed $`D^+`$ or $`D^0`$ candidate. The selection of muons and electrons is described in detail in . For this analysis, lepton candidates are required to have a momentum of at least 2.0 GeV/$`c`$ for electrons and 2.5 GeV/$`c`$ for muons.

The $`D^+`$ and $`D^0`$ candidates are reconstructed from charged tracks and $`\pi ^0`$’s that form an angle of less than $`45^{}`$ with the lepton candidate. These charged tracks are also required to intersect a cylinder of radius 2 cm and half-length 4 cm centered on the nominal interaction point, to have at least 4 hits in the TPC, a polar angle $`\theta `$ such that $`|\mathrm{cos}\theta |<0.95`$ and a transverse momentum greater than 200 MeV/$`c`$.

Photons and $`\pi ^0`$’s are identified in the ECAL. The four-momenta of $`\pi ^0`$ candidates are computed by adding the photon momenta when the $`\gamma \gamma `$ invariant mass is consistent with the $`\pi ^0`$ mass. The energy of the $`\pi ^0`$ is then recomputed using the kinematical constraint of the $`\pi ^0`$ mass . The energy resolution achieved is about 6.5%, almost independent of the energy. For reconstructing the $`D^0`$ candidates, only $`\pi ^0`$’s with momenta greater than 2 GeV/$`c`$ are used.

For all charged kaons used to reconstruct the $`D^0`$ it is required that the $`dE/dx`$ be within $`3\sigma `$ of that expected from a kaon, except for the $`D^+D^0\pi ^+`$, $`D^0K^{}\pi ^+`$decay channel. The kaon has the same sign as the lepton coming from the $`B`$ semileptonic decay, therefore this charge correlation is required.

$`K_\mathrm{S}^0`$ candidates are reconstructed from pairs of oppositely charged tracks. The two tracks are required to be inconsistent with originating from the interaction point and the $`K_\mathrm{S}^0`$ candidate is rejected if the measured mass is more than $`2\sigma `$ ($`\pm 10`$MeV/$`c^2`$) from the nominal $`K_\mathrm{S}^0`$ mass.

Tracks coming from the $`D^0`$ decays are required to form a common vertex with a $`\chi ^2`$ confidence level greater than 1%. In the $`D^0`$ decay channels that do not contain a $`\pi ^0`$ in the final state, if more than one combination satisfies the selection criteria, the one with the smallest $`\chi ^2`$ of the $`D`$ vertex fit is selected. For the decay modes where the multiple candidates originate from different $`\pi ^0`$ combinations, for a given detected lepton, the candidate with $`m_{D^0\pi }m_{D^0}`$ closest to the nominal value is chosen in case of $`D^{}`$ selection, otherwise the highest momentum $`D^0`$ candidate is chosen.

In order to reject the background from charm, to improve the signal to background ratio and to ensure well-measured decay lengths, additional selection criteria are applied to all the subsamples. The invariant mass of the $`D^{()}\mathrm{}`$ system is required to be greater than 3 GeV/$`c^2`$, where $`D^{()}`$ can be $`D^+`$ or $`D^0`$. This cut significantly reduces the charm background while keeping $`85\%`$ of the signal. To exploit the high precision of the silicon vertex detector, the lepton track and at least two tracks from the $`D^0`$ decay are required to have at least one VDET hit in both the $`r\varphi `$ and $`z`$ projections.

### 3.1 The $`𝑫^{\mathbf{}}`$-lepton selection

$`D^+`$ candidates are identified via the decay $`D^+D^0\pi ^+`$, followed by $`D^0K^{}\pi ^+`$, $`D^0K^{}\pi ^+\pi ^{}\pi ^+`$, $`D^0K^{}\pi ^+\pi ^0`$ or $`D^0K_\mathrm{S}^0\pi ^+\pi ^{}`$. The difference in mass between the $`D^+`$ and $`D^0`$ candidates must lie within 1.5 MeV/$`c^2`$ (approximately two standard deviations of the experimental resolution) of the nominal value of 145.4 MeV/$`c^2`$.

In the $`D^0K^{}\pi ^+`$ channel, the $`D^0`$ momentum $`p_{D^0}`$ is required to be greater than 5 GeV/$`c`$.

For the $`D^0K^{}\pi ^+\pi ^{}\pi ^+`$ channel $`p_{D^0}`$ must be greater than 8 GeV/$`c`$, and at least two of the $`D^0`$ decay tracks must have $`p>1`$ GeV/$`c`$.

In the case of the $`D^0K^{}\pi ^+\pi ^0`$ channel, the momentum of the reconstructed $`D^0`$ is required to be greater than 10 GeV/$`c`$, and the two charged tracks in the decay are required to have $`p>`$ 0.5 GeV/$`c`$. Futhermore the decay kinematics are required to be consistent with one of the three resonant decays: $`D^0K^{}\rho ^+`$$`D^0K^{}\pi ^+`$, $`D^0\overline{K}^0\pi ^0`$. For each decay hypothesis the mass of the resonant particle and the helicity angle $`\theta _H`$ are calculated.111 The helicity angle $`\theta _H`$ is defined as the angle between the scalar particle and one of the decay products of the vector particle, calculated in the rest frame of the vector particle.  If the mass is consistent with the nominal resonance mass value, within its natural width, and if $`|\mathrm{cos}\theta _H|>0.4`$, the candidate is considered to be consistent with the resonant decay hypothesis. For these resonant decays a $`\mathrm{cos}^2\theta _H`$ distribution is expected.

Finally, in the $`D^0K_\mathrm{S}^0\pi ^+\pi ^{}`$ channel, the momentum of the reconstructed $`D^0`$ is required to be greater than 4 GeV/$`c`$. The same technique as in the $`D^0K^{}\pi ^+\pi ^0`$ channel is used to tag the $`K^{}`$ resonance. The sign of the resonance is used to distinguish between $`D^0`$ and $`\overline{D}^0`$.

### 3.2 The $`𝑫^\mathrm{𝟎}`$-lepton selection

The $`D^0\mathrm{}^{}`$ sample consists of events with a lepton and a $`D^0`$ candidate, where the $`D^0`$ is not the decay product of a reconstructed $`D^+`$. The $`D^0`$ candidates are identified using the same decay modes as for the $`D^+\mathrm{}^{}`$sample. For this sample, the background is larger because the $`D^+`$$`D^0`$ mass difference criterion is not applicable.

For all the decay modes selected in the $`D^0\mathrm{}^{}`$sample, a search for the additional pion is performed to reject $`D^0`$ candidates coming from $`D^+D^0\pi ^+`$. If a pion candidate is found yielding a $`D^+`$$`D^0`$ mass difference within 6 MeV/$`c^2`$ of the nominal value, the $`D^0\mathrm{}^{}`$ candidate is rejected. The efficiency for reconstructing the additional pion and rejecting $`D^0`$’s coming from $`D^+`$ decays is found to be 86% from the Monte Carlo simulation.

In the $`D^0K^{}\pi ^+`$ channel, the reconstructed $`D^0`$’s are required to have $`\text{p}_{D^0}>8\text{GeV/}c`$, $`p_K>2`$ GeV/$`c`$ and $`p_\pi >1.5`$ GeV/$`c`$. In the $`D^0K^{}\pi ^+\pi ^{}\pi ^+`$ channel the $`D^0`$ momentum is required to exceed $`12`$ GeV/$`c`$ while the kaon and the three pions must have momenta greater than 2 and 1 GeV/$`c`$, respectively.

The cuts on the kaon and pion momenta are tightened for the $`D^0K^{}\pi ^+\pi ^0`$ channel, to $`p_K>3`$ GeV/$`c`$ and $`p_\pi >2`$ GeV/$`c`$, while the reconstructed $`D^0`$ must have $`\text{p}_{D^0}>12\text{GeV/}c.`$ The same cuts on the three resonances are used as in the $`D^+\mathrm{}^{}`$event sample.

Finally, for the $`D^0K_\mathrm{S}^0\pi ^+\pi ^{}`$ channel a cut is applied of at least $`1.5`$ GeV/$`c`$ on the kaon momentum and $`1.0`$ GeV/$`c`$ on the two pion momenta. The momentum of the reconstructed $`D^0`$ candidates is required to be greater than 10 GeV/$`c`$. The same technique is used to tag the $`K^{}`$ resonance as in the $`D^+\mathrm{}^{}`$ event sample.

### 3.3 $`𝑩`$ meson reconstruction

The $`B`$ decay vertex position is estimated by vertexing the reconstructed $`D^0`$ track with the lepton. Events are rejected if the $`B`$ vertex fit gives a $`\chi ^2`$ probability less than 1%.

The $`D^0`$ candidate mass spectra for the four subsamples in the $`D^+\mathrm{}^{}`$ event selection are shown in Fig. 1, and in the $`D^0\mathrm{}^{}`$ event selection in Fig. 2. The fitted curves consist of a Gaussian function for the signal plus a linear background. For the $`D^0K^{}\pi ^+\pi ^0`$ in the $`D^0\mathrm{}^{}`$sample a Gaussian tail is used to describe the background; masses below 1.65 GeV/$`c^2`$ are excluded to avoid the broad enhancement due to the missing $`\pi ^0`$ in the $`D^0K^{}\pi ^+\pi ^0\pi ^0`$ decay. The fitted $`D^0`$ width and the fitted number of signal and background events within a window of $`\pm 2\sigma `$ around the fitted mass are shown in Table 1 for the different samples.

Events reconstructed within two standard deviations of the fitted $`D^0`$ mass are selected for the lifetime analysis, resulting in 1880 $`D^+\mathrm{}^{}`$ and 2856 $`D^0\mathrm{}^{}`$ candidates. The decay length is calculated for these events by reconstructing the primary and $`B`$ decay vertices in three dimensions. The primary vertex reconstruction algorithm  applied to simulated $`b\overline{b}`$ events yields an average resolution of $`50\mu \text{m}\times 10\mu \text{m}\times 60\mu \text{m}`$ (horizontal $`\times `$ vertical $`\times `$ beam direction).

The distance between the primary and $`B`$ decay vertices is projected onto the direction defined by the momentum of the $`D^{()}\mathrm{}`$ system. The uncertainty on the flight direction due to the missing neutrino induces a negligible error on the decay length. The resolution on the $`B`$ decay length is on average 250 $`\mu `$m, compared with an average $`B`$ decay length of about $`2.5`$ mm.

For the $`D^0K^{}\pi ^+\pi ^0`$ channel, the $`\pi ^0`$ momentum is included when extrapolating the neutral $`D^0`$ track backwards to form the $`B`$ vertex. In the case of $`D^+\mathrm{}^{}`$ events, the soft pion from the $`D^+`$ decay does not improve the resolution on the $`B`$ decay length and is therefore not used in the reconstruction of the $`B`$ vertex.

Because the selected decays contain an undetected neutrino, the $`B`$ momentum is not known precisely and is reconstructed using an energy flow technique as described in . A further correction is applied by evaluating for Monte Carlo events the $`\kappa `$ distribution, defined as

$$\kappa =\frac{(\beta \gamma )_{\mathrm{reco}}}{(\beta \gamma )_B}.$$
(1)
Because this distribution depends on the selection criteria applied, separate $`\kappa `$ distributions are calculated for each subsample. The momentum resolution obtained is between 8 and 10%, depending on the decay channel.

## 4 Lifetime measurement

An unbinned likelihood fit for the lifetimes is performed. For each event the probability of observing a proper time $`t`$ given the lifetime $`\tau `$ is calculated: the probability density function $`F(t,\sigma _t,\tau )`$ is obtained by convoluting an exponential distribution with the properly normalised $`\kappa `$ distribution and with a Gaussian function which takes into account the resolution on the decay length.

Both the $`D^+\mathrm{}^{}`$ and $`D^0\mathrm{}^{}`$ samples contain a mixture of $`\overline{B}^0`$ and $`B^{}`$ decays and the $`B^{}`$/$`\overline{B}^0`$ mixture in the samples depends on the ratio of the lifetimes, as discussed in Section 4.2. Therefore, to measure the $`\overline{B}^0`$ and $`B^{}`$ lifetimes a simultaneous maximum likelihood fit is performed to all the events. The likelihood function contains three components for each sample and is written as

$``$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N_D^{}\mathrm{}}{}}}\left[f_{}^{}(\tau _B^{}/\tau _{\overline{B}^0})F(t_i,\sigma _i,\tau _B^{})+f_0^{}(\tau _B^{}/\tau _{\overline{B}^0})F(t_i,\sigma _i,\tau _{\overline{B}^0})+f_{\mathrm{bkg}}^{}F_{\mathrm{bkg}}^{}(t_i)\right]`$ (2)
$`\times `$ $`{\displaystyle \underset{i=1}{\overset{N_{D^0\mathrm{}}}{}}}\left[f_{}^0(\tau _B^{}/\tau _{\overline{B}^0})F(t_i,\sigma _i,\tau _B^{})+f_0^0(\tau _B^{}/\tau _{\overline{B}^0})F(t_i,\sigma _i,\tau _{\overline{B}^0})+f_{\mathrm{bkg}}^0F_{\mathrm{bkg}}^0(t_i)\right].`$
The coefficients $`f_{}^{}`$ and $`f_0^{}`$ are the fractions of the $`D^+\mathrm{}^{}`$ sample arising from $`B^{}`$ and $`\overline{B}^0`$ decays, respectively. Similarly, $`f_{}^0`$ and $`f_0^0`$ are the fractions of $`B^{}`$ and $`\overline{B}^0`$ decays in the $`D^0\mathrm{}^{}`$ sample. The coefficients $`f_{\mathrm{bkg}}^{}`$ and $`f_{\mathrm{bkg}}^0`$ are the background fractions of the samples, while the functions $`F_{\mathrm{bkg}}^{}(t)`$ and $`F_{\mathrm{bkg}}^0(t)`$ are their normalised proper time distributions.

### 4.1 Backgrounds

Background contamination arises from the following sources:

1. combinatorial background, i.e. candidates with a fake $`D^{()}`$;
2. the process $`\overline{B}D_s^{}D^{()}X`$, followed by $`D_s^{}\mathrm{}^{}X`$, giving rise to a real $`D^{()}`$ and a real lepton;
3. a real $`D^{()}`$ meson accompanied by a fake or nonprompt lepton, from $`Zb\overline{b}`$ or $`Zc\overline{c}`$ events.

Source (1) is the dominant background. Its contribution is determined from a fit to the $`D^0`$ mass distributions, and its magnitude is given in Table 1 for the various subsamples. The proper time distribution for this source is determined from the data by selecting events from the high mass sideband of the $`D^0`$ peak. The same selection criteria described in Section 3 are applied to the background samples, except that the requirement on the $`D^+`$$`D^0`$ mass difference in the case of the $`D^+\mathrm{}^{}`$ events is removed to increase the statistics. A function consisting of a Gaussian plus positive and negative exponential tails is used to describe these data.

The contribution from source (2) is calculated from the measured branching ratios for this process  and a Monte Carlo simulation to determine the detection efficiency; it accounts for a contamination which is about 2–3% of the sample, depending on the channel.

The background from source (3) is estimated from the measured hadron-lepton misidentification probabilities  and the measured inclusive $`D^0`$ and $`D^+`$ rates. An independent estimate is obtained using wrong-sign ($`D^+\mathrm{}^+`$ or $`D^0\mathrm{}^+`$) events, and is found to be consistent. This background source contributes between $`2`$ and $`5\%`$ of the sample, depending on the channel. To characterize the proper time distribution of this background, it is further subdivided into three distinct components: fake leptons coming from the primary vertex in $`c\overline{c}`$ and $`b\overline{b}`$ events, and fake leptons coming from a decaying $`b`$ hadron.

The proper time distribution for sources (2) and (3) are determined from simulated events.

### 4.2 Sample compositions

Both the $`D^+\mathrm{}^{}`$ and $`D^0\mathrm{}^{}`$ samples contain a mixture of $`\overline{B}^0`$ and $`B^{}`$ decays. In order to estimate the cross contamination, the individual semileptonic branching ratios of the $`\overline{B}^0`$ and $`B^{}`$ must be determined. The evaluation follows the same procedure as given in the appendix of . An important input in this evaluation is the knowledge of the branching ratios for the decay modes $`BD^{}\pi \mathrm{}\nu `$, where $`D^{}`$ represents any of the $`\text{p}\mathrm{wave}`$ resonances as well as the nonresonant $`D^{()}\pi `$ states. The most recent ALEPH and DELPHI results  for both the resonant and the nonresonant components are used in the calculation, leading to a significant reduction in the resulting uncertainty compared to the previous analysis of . The $`\overline{B}^0`$ and $`B^{}`$ content of the two samples are calculated using as input the measured values of the branching ratios given in Table 2. The sample composition is then calculated by considering the $`\overline{B}^0`$ and $`B^{}`$ decay channels that contribute to the $`D^+\mathrm{}^{}`$ and $`D^0\mathrm{}^{}`$ samples , taking into account the probability of $`0.147\pm 0.015`$ that a $`D^+\mathrm{}^{}`$ event is mistakenly reconstructed as a $`D^0\mathrm{}^{}`$ event.

As a consequence of this procedure, the coefficients $`f_{}^{}`$, $`f_0^{}`$, $`f_{}^0`$ and $`f_0^0`$ appearing in the likelihood function (Eq. 2) depend on the lifetime ratio. For equal lifetimes, about $`85\%`$ of the $`B`$ decays in the $`D^+\mathrm{}^{}`$ sample are attributed to $`\overline{B}^0`$, while about $`80\%`$ of the $`D^0\mathrm{}^{}`$ sample $`B`$ decays come from $`B^{}`$.

### 4.3 Fit results

The fit to the proper time distributions of the $`D^+\mathrm{}^{}`$ and $`D^0\mathrm{}^{}`$ events is performed to determine the two free parameters $`\tau _{\overline{B}^0}`$ and $`\tau _B^{}`$. The values obtained are

$`\tau _{\overline{B}^0}`$ $`=`$ $`1.518\pm 0.053\mathrm{ps},`$
$`\tau _B^{}`$ $`=`$ $`1.648\pm 0.049\mathrm{ps},`$
where the errors are statistical only. The statistical correlation coefficient is $`0.35`$. The ratio of the lifetimes is

$$\tau _B^{}/\tau _{\overline{B}^0}=1.085\pm 0.059,$$

taking into account the correlation.

The proper time distributions for the two samples are shown in Fig. 3, with the results of the fit overlaid.

### 4.4 Systematic uncertainties

The sources of systematic uncertainty are discussed in the following, and the estimated errors are summarized in Table 3.

The uncertainty in the $`B`$ momentum reconstruction is dominated by the uncertainty in the $`\kappa `$ distribution. The effects that can modify the $`\kappa `$ function with a consequence on the extracted lifetimes are studied. The first is due to the $`D^{()}\pi `$ content of semileptonic $`B`$ decays, which affects the $`\kappa `$ function due to the smaller phase space available and the softer momentum spectrum of the reconstructed $`D^{()}\mathrm{}`$ system. The fractions of $`D^{()}\pi `$ decays are varied in the simulation within experimental errors. Other effects studied are the dependence of the $`\kappa `$ distribution on the analysis cuts and the $`b`$ fragmentation function. These two effects combined give a relative uncertainty of about $`2\%`$ in the $`B`$ momentum determination. Uncertainties in the $`\kappa `$ distribution propagate in the same way to $`\tau _B^{}`$ and $`\tau _{\overline{B}^0}`$, and therefore have a small effect in the ratio.

Uncertainties in the background fractions and proper time distributions are considered. The combinatorial background, estimated in the data from fits to the $`D^0`$ mass spectra, is reported in Table 1 for each individual channel. The systematic uncertainty due to this source is estimated by varying the combinatorial background in the fit for each channel within its statistical errors in turn and taking the sum in quadrature. The contamination from physics background is evaluated using the simulation, by varying the fraction within the estimated uncertainty and repeating the fits. The uncertainty on the total background level is calculated by combining the two above sources.

The parameters describing the background proper time distributions are varied within their uncertainties. Background proper time distributions are parametrized using different methods, to check for possible systematic bias. Different background samples are selected by varying the sideband regions, adding events from the lower sideband, which are excluded in the lifetime determination, or using events with wrong-sign correlations. Some cuts in the selection are varied to check the stability of the parametrizations. The shapes extracted from real data are compared with those extracted from Monte Carlo events. The resulting differences in the fitted lifetimes are used to estimate the systematic uncertainty.

For the remaining systematic errors, which are small compared to the statistical ones, the correlation between the lifetimes is not propagated into the ratio.

The systematic uncertainty due to the sample compositions is determined by varying the branching fractions of Table 2 within $`\pm 1\sigma `$ from the central values. In addition the uncertainty due to the assumption of isospin conservation has been estimated, allowing for a variation of 20% relative to the exact symmetry.

As explained in detail in , a small fraction of four-body decays where an extra pion is produced in the decay of the B meson contributes to the sample composition, and is characterized by a different selection efficiency with respect to the real signal. The relative efficiencies

$$\frac{ϵ(BD^{()}\pi \mathrm{}\nu )}{ϵ(BD^{()}\mathrm{}\nu )}=0.75\pm 0.10(\text{for }\overline{B}^0)0.64\pm 0.10(\text{for }B^{})$$

therefore enter into the calculation and this uncertainty (coming from the limited Monte Carlo statistics) is propagated to the measured lifetimes.

The parameters of the decay length resolution function are varied within their errors to estimate the resulting uncertainty .

## 5 Conclusions

The lifetimes of the charged and neutral $`B`$ mesons have been measured with the full statistics collected by the Aleph detector at and around the Z peak energy. The data sample was recently reprocessed, achieving improved tracking performance. A maximum likelihood fit to the proper time distributions of 1880 $`D^+\mathrm{}^{}`$ and 2856 $`D^0\mathrm{}^{}`$ candidates yields the following results for the $`\overline{B}^0`$ and $`B^{}`$ lifetimes and their ratio:

$`\tau _{\overline{B}^0}`$ $`=`$ $`1.518\pm 0.053\pm 0.034\mathrm{ps},`$
$`\tau _B^{}`$ $`=`$ $`1.648\pm 0.049\pm 0.035\mathrm{ps},`$
$`\tau _B^{}/\tau _{\overline{B}^0}`$ $`=`$ $`1.085\pm 0.059\pm 0.018,`$
where the first error is statistical and the second is systematic.

These results supersede the corresponding ones of . Averaging with the results of the other methods presented in  the combined values are:

$`\tau _{\overline{B}^0}`$ $`=`$ $`1.496\pm 0.048\pm 0.033\mathrm{ps},`$
$`\tau _B^{}`$ $`=`$ $`1.644\pm 0.048\pm 0.034\mathrm{ps},`$
$`\tau _B^{}/\tau _{\overline{B}^0}`$ $`=`$ $`1.104\pm 0.057\pm 0.019.`$
## Acknowledgments

We thank our colleagues in the accelerator divisions for the excellent performance of LEP. Thanks also to the many engineering and technical personnel at CERN and at the home institutes for their contributions to the performance of the Aleph detector. Those of us from non-member states thank CERN for its hospitality.