text:
# 1 Introduction

## 1 Introduction

An outstanding problem in high energy nuclear scattering is whether the hot and dense matter formed equilibrates to briefly form a plasma of deconfined quarks and gluons–the quark gluon plasma. This problem is of great topical interest, with collisions already taking place at the Relativistic Heavy Ion Collider (RHIC) and scheduled to take place several years hence at the Large Hadron Collider (LHC).

Whether or not a quark gluon plasma is formed depends strongly on the highly non–equilibrium initial distributions of partons formed immediately after the collision. Clearly, these distributions must influence the subsequent evolution of the system towards equilibrium. Furthermore, the problem is complicated by the rapid expansion of the system as a whole. The magnitude of the collision induced relaxation time relative to the expansion time is what determines whether equilibrium is indeed reached.

The study of equilibration in relativistic heavy ion collisions is as old as the subject itself. However, very few ab initio studies exist that attempt to follow the evolution of the system, all the way, from the first instants of the collision to equilibrium . This was so because it was not known how to treat, in a self–consistent manner, the small Bjorken $`x`$ “wee” parton modes that are responsible for particle production at central rapidities. These small $`x`$ modes provide the initial conditions for the space–time evolution of the partonic matter formed in heavy ion collisions.

Because the occupation number of small $`x`$ modes in the nuclear wavefunctions is large, it was shown that classical methods could be used to compute their distributions . This classical effective theory is, in exact analogy to a spin glass, a color glass condensate , and is characterized by a bulk scale $`Q_s`$ – the momentum scale at which gluon distributions saturate. Analytical expressions for the parton distributions were obtained in Refs.  and . Subsequently, quantum corrections to the classical “non-Abelian Weizsäcker–Williams” fields were computed , and renormalization group methods  devised to study how the classical parton distributions in the nuclei changed with energy (or equivalently, with $`x`$). The scale $`Q_s`$ now depends on $`x`$. This dependence is represented by a line in the $`x`$$`Q^2`$ plane – it separates the saturated, non–linear regime of QCD at high parton densities from that of linear QCD evolution. That $`Q_s`$ is a function of $`x`$, and grows as one goes to smaller $`x`$, or equivalently, to higher energies, will only be implicit in this work.

In the classical effective field theory approach, the problem of initial conditions can be formulated as the problem of finding solutions of the Yang–Mills equation with initial conditions given by the classical fields of each of the nuclei before the collision . Since analytical expressions exist for the classical fields of the nuclei before the collision, the initial conditions are fully determined.

Perturbative solutions of the Yang–Mills equations, which describe classical gluon production to lowest order, have been discussed by several authors . These were found to be infra–red divergent. Within the Yang–Mills approach, a fully non–perturbative treatment is therefore necessary. Non–perturbative, numerical solutions of the Yang–Mills equations have been found recently . In particular, one now knows the initial number and energy distribution of gluons after a collision . An idealized form of the initial gluon distribution, in terms of the saturation scale of the color glass condensate $`Q_s`$, was given recently by Mueller . These idealized distributions are sufficient for the purposes of this paper. We will reserve a more quantitative analysis using the initial distributions of Ref.  for a future work .

The initial partonic system is completely out of equilibrium. The subsequent scattering and evolution of the system towards equilibrium was studied by Mueller . In Ref. , he showed that, under the assumption that the system is undergoing boost invariant expansion, the evolution of single particle distributions could be described by a non–linear Landau equation. At very early times, this equation can be linearized ***to a Fokker-Planck equation ., and studied analytically. However, the analytical approximations soon break down, and the evolution of the system cannot be followed analytically all the way to equilibrium.

In this paper, we numerically solve the Landau equation proposed by Mueller. (Along the way, we compare our results to Mueller’s analytical results and show quantitatively where the analytical approach breaks down.) We are thus able to follow the evolution of the system all the way to equilibrium. We study the dependence of the equilibration time, the initial temperature, and the chemical potential on the saturation scale $`Q_s`$ of the nuclear wavefunction, and on $`\alpha _S`$. (Note: one may estimate $`Q_s1`$ GeV at RHIC and $`Q_s23`$ GeV at LHC .) We discuss the dynamical screening of the collinear divergence arising from small angle scattering. Our results for equilibration are obtained, primarily, in a particular model of dynamical screening .

To conclude whether local kinetic equilibrium is attained in a system undergoing boost invariant 1+1–dimensional expansion, one should compare the equilibration time $`t_{equil}`$ to the hydrodynamic expansion time $`t_{hydro}`$. Typically, $`t_{hydro}R/c_0`$, where $`R`$ is the radius of the nucleus, and $`c_0`$ is the speed of sound in the system. For an ideal ultrarelativistic gas, $`c_0=1/\sqrt{3}`$. Our qualitative results suggest that for the energies of interest at RHIC and LHC, $`t_{equil}t_{hydro}`$, thereby indicating that favorable conditions may exist for the formation of a hot gluon plasma at RHIC and at LHC.

In studying the likelihood of equilibration, we have only considered number conserving $`22`$ processes. Naively, $`23`$ processes will be suppressed by a power of $`\alpha _S`$, and may be considered sub–leading. However, $`23`$ processes may, in principle, be much more efficient than $`22`$ processes in driving the system towards thermal equilibrium . Whether this is indeed the case is a dynamical question we will not address in this work. Regardless, $`22`$ processes should set an upper bound on the equilibration time. Furthermore, it is likely that the numerical techniques developed here can be adapted to quantitatively study the effects of number changing processes.

The paper is organized as follows. We begin in section 2 by discussing the initial conditions for gluon production in nuclear collisions. After the collision, particle production is described by the space–time evolution of classical gauge fields. At late times, when the system is dilute, one can define partons, and study their interactions in a kinetic approach. The kinetic transport equation derived by Mueller is described in section 3. In the following section, we discuss the problem of dynamical screening of collinear divergences, and describe the particular model of screening that is employed in this work. In section 5, we first briefly outline the numerical method, and compare results for a linearized (Fokker–Planck) version of the Landau equation with Mueller’s analytical results. We next discuss numerical results from the solution of the full Landau equation. Various measures of equilibration are described and computed as a function of the saturation scale $`Q_s`$ and the coupling $`\alpha _S`$. We summarize our results in the final section and comment on their experimental ramifications. We will also discuss the various uncertainties that need to be further quantified, and their likely impact on equilibration in high energy heavy ion collisions. An appendix contains a detailed description of the numerical method employed to solve the Landau equation.

## 2 The initial conditions for parton evolution

The problem of initial conditions in heavy ion collisions is as follows. Before the collision, the nucleus, as viewed in a particular gauge, is a coherent superposition of various Fock modes containing differing numbers of partons For an illuminating exposition of this point of view, see Ref. .. Partons in Fock modes containing a large number of partons, $`|qqq\mathrm{}qgggg\mathrm{}q\overline{q}`$, have small values of $`x`$–in the language of the parton model, the momentum of the nucleus is shared among a large number of constituents. Each parton has a $`+`$ momentum $`k^+`$, and $`x=k^+/P^+`$, where $`P^+`$ is the momentum of the nucleus. The large multiplicities at central rapidities–which correspond to small $`x`$–are obtained when these small $`x`$ Fock modes “go on-shell”. One should note that even though the small $`x`$ modes carry small longitudinal momenta, their transverse momenta $`p_t`$ can be large.

These small $`x`$ partons are highly de–localized in the longitudinal $`x^{}`$ direction if $`P^+`$ is the momentum of the nucleus – they have large wavelengths relative to the Lorentz contracted width of the nucleus. In addition, these virtual fluctuations are very short–lived. Any attempt to treat them in a naively classical transport approach is therefore problematic. An approach based on pQCD collinear factorization is also inadequate unless the partons also have very large transverse momenta. This is because such a picture is predicated on convolving the probabilities of parton distributions in the two nuclei with the elementary parton–parton scattering cross-section. It cannot describe coherence effects which are important at small $`x`$. (For large $`p_t`$, these effects are suppressed.)

It was realized some time ago that the small $`x`$ Fock states in nuclei responsible for multi–particle production at central rapidities are states of high occupation number. As a consequence, their distributions in the nuclei are described by a classical effective field theory (EFT) . The classical distributions in a single nucleus can be solved and an analytical form for the distributions obtained . The classical gluon distribution falls off as $`1/k_t^2`$ at large transverse momentum $`k_t`$ but saturates at smaller $`k_t`$–growing only slowly as $`\mathrm{ln}(\alpha _S\mu /k_t)/\alpha _S`$. The infrared structure of the EFT is analogous to that of a spin glass in condensed matter systems –it is thus a colored glass condensate .

The EFT approach was used in Ref.  to treat the problem of initial conditions in nuclear collisions. It consists of solving the Yang–Mills equations in the forward light cone, after the collision, with initial conditions given by the (known) gauge fields of the nuclei before the collision.

Gluons at central rapidities are produced classically. (Boost invariance is assumed.) Perturbative computations of the transverse momentum distributions of classically produced gluons (lowest order in $`\alpha _S`$ and lowest non–trivial order in $`\alpha _SQ_s/p_t`$) were first computed in Ref. , and later by several authors . In these computations, the number distribution is found to be infra-red divergent. They agree with the lowest order mini–jet predictions at large $`p_t`$ §§§For detailed comparisons, see the papers by Gyulassy and McLerran, and by Guo, in Ref. .

When $`p_tQ_s`$, all orders in $`Q_s/p_t`$ contribute equally. The perturbative computations, to lowest order in $`Q_s/p_t`$ are therefore not sufficient. Thus, even at lowest order in $`\alpha _S`$, one needs to resum all orders in $`Q_s/p_t`$. An analytic solution to this non–perturbative problem has not been found. However, the problem was recently formulated as a 2+1–dimensional classical effective field theory and solved numerically (for an SU(2) gauge theory) on a 2–dimensional transverse lattice . The evolution of gauge fields is computed in real time, and at late times, the energy  and the number  of produced gluons at central rapidities can be computed. The initial distribution of gluons, from our numerical simulations , is as follows. At low values of $`p_t`$, $`p_tQ_s`$ it is remarkably like the Bose–Einstein distribution for a gas of massive particles in two dimensions. The gluons acquire a screening mass due to strong non–linear interactions – this also renders the distributions infrared finite. At large $`p_t`$, $`p_tQ_s`$, the distribution is a power law, $`Q_s^4/p_t^4`$. The Born–pQCD prediction for the gluon distribution is therefore obtained, as expected, in the large $`p_t`$ limit.

The gluons produced at central rapidities are completely out of equilibrium. The study of the evolution of this system towards equilibrium was initiated by Mueller . For simplicity, Mueller chose idealized initial conditions–a constant instead of a Bose–Einstein distribution for $`p_tQ_s`$, and zero for $`p_tQ_s`$–namely, the theta function $`\theta (Q_s^2p_t^2)`$. In this work, we will consider only these idealized distributions of Mueller. In a later work , we will consider more realistic distributions computed very recently from the lattice simulations of Ref. . In general, the initial spatial distribution per unit rapidity is

$`{\displaystyle \frac{1}{L^2}}{\displaystyle \frac{dN}{d\eta }}=c{\displaystyle \frac{N_c^21}{4\pi ^2\alpha _SN_c}}Q_s^2,`$ (1)
where $`L^2`$ is the transverse area, $`\eta `$ is the space–time rapidity, and $`c`$, determined non–perturbatively, is a weak function of $`Q_sR`$. For an SU(2) gauge theory, it was computed recently in Ref.  to be $`c=1.29\pm 0.09`$ in the regime of interest.

The initial distributions are highly anisotropic, with the produced partons having zero longitudinal momentum $`p_z`$. Assuming boost invariance, and using the relation

$`f(x,p)={\displaystyle \frac{(2\pi )^3}{2(N_c^21)}}{\displaystyle \frac{dN}{d^3pd^3x}},`$ (2)
for the single particle distribution $`f(x,p)`$, one finds at late times that The $`\delta `$–function distribution is of order $`1/\mathrm{\Delta }p_z`$. We shall see later, from the early time analytical solution of the Fokker–Planck equation, that the longitudinal momentum changes little in the time in takes for the occupation number to decrease below unity.

$`f(x,p)={\displaystyle \frac{c}{\alpha _SN_c}}{\displaystyle \frac{1}{t}}\delta (p_z)\theta (Q_s^2p_t^2).`$ (3)
At late times, the occupation number of partons in the transverse plane becomes small, and their evolution can no longer be followed on the lattice. However, when the occupation number is not too large, transport theory can be applied to study the further evolution of the system. The initial condition for this evolution is given by the single particle distribution in Eq. 3.

Using the results of Ref. , which determined the time at which gluons come on-shell after a collision, we can estimate the time at which small angle scattering between the gluons can be described by the Boltzmann equation. The formation times, using $`Q_s=1`$ GeV for RHIC and $`Q_s=2`$$`3`$ GeV for LHC, are

$$t_i1.40\mathrm{GeV}^1\mathrm{for}\mathrm{RHIC},t_i0.62\mathrm{GeV}^1\mathrm{for}\mathrm{LHC}$$
(4)
Next, we suppose that a transport theory based treatment of gluon scattering is applicable when the gluon occupation number is less than unity. Letting $`d\eta `$ be $`d\eta =dz/t_i=1`$ or $`dz=t_i`$, and using $`\frac{dN}{d^3xd^3k}`$ from  at $`z=0`$, we find

$$t_0=\frac{c}{\alpha _SN_c}t_i.$$
(5)
Here $`t_0>t_i`$ is the time where the gluon occupation number is dilute enough that their subsequent interactions can be described as the small angle scattering of on-shell gluons.

## 3 The Landau transport equation

At times $`t>t_0`$, where $`t_0`$ was defined in Eq. 5, a Boltzmann–like transport equation is appropriate for describing the late–time evolution of the highly anisotropic initial gluon distribution discussed in the previous section. We will show below that, for central high energy heavy ion collisions, assuming boost invariance, the transport equation reduces to a Landau–type transport equation .

In the problem of interest, the typical scale of spatial variations is large compared to the typical scale of gluon–gluon scattering. The evolution of the system can therefore be described by a local Boltzmann–like kinetic equation for the single particle gluon distribution $`f(x,p,t)`$, which gives the density of excitations of momentum $`\stackrel{}{p}`$ at a point $`(\stackrel{}{x},t)`$. The Boltzmann equation is

$`{\displaystyle \frac{f}{t}}+\stackrel{}{v}\stackrel{}{}f=C(f),`$ (6)
where $`C(f)`$ is the collision integral which represents the change, due to collisions, per unit time, in the number of particles per unit phase space volume.

Let us first consider the left hand side of the above equation. We will assume that the transverse dimensions of the collision volume in central heavy ion collisions are sufficiently large that the initial expansion of the system is one–dimensional. The distribution $`f`$ then only depends on $`z`$, the co–ordinate corresponding to the collision axis. Since it is independent of the transverse co–ordinate $`x_t`$, Eq. 6 has the form

$`{\displaystyle \frac{f(\stackrel{}{p},z,t)}{t}}+v_{p_z}{\displaystyle \frac{f(\stackrel{}{p},z,t)}{dz}}=C\left[f(\stackrel{}{p},z,t)\right],`$ (7)
where $`v_{p_z}=p_z/|\stackrel{}{p}|`$. A further simplification follows from the assumption that the central rapidity region is approximately Lorentz invariant under boosts . As was pointed out by Baym many years ago , this assumption greatly simplifies the problem of solving Eq. 7 since it relates the distribution function at different $`z`$’s in the central region. It therefore suffices to compute $`f`$ in the zero–rapidity slice alone. Because $`f`$ is a scalar under Lorentz transformations, it satisfies the relation

$`f(p_t,p_z,z,t)=f(p_t,p_z^{},t).`$ (8)
Here $`p_z^{}=\gamma (p_zup)`$, the transformation velocity $`u=z/t`$, and $`\gamma =(1u)^{1/2}=t/\tau `$, where the proper time $`\tau =\sqrt{t^2z^2}`$.

Computing $`f/z`$ using Eq. 7, the Lorentz transformation relation, $`\tau /z|_{z=0}=0`$, and $`p_z^{}/z|_{z=0}=p/t`$, one finds  at $`z=0`$,

$`v_{p_z}{\displaystyle \frac{f}{z}}={\displaystyle \frac{p_z}{t}}{\displaystyle \frac{f}{p_z}}.`$ (9)
The Boltzmann equation now reduces to

$`{\displaystyle \frac{f(p_t,p_z,t)}{t}}|_{p_zt}\left({\displaystyle \frac{}{t}}{\displaystyle \frac{p_z}{t}}{\displaystyle \frac{}{p_z}}\right)f(p_t,p_z,t)=C\left[f(p_t,p_z,t)\right].`$ (10)
Note that, as a consequence of our manipulations, assuming uniform distributions in the transverse direction and boost invariance in the longitudinal direction, the single particle distribution $`f(x,p,t)`$ is now expressed as a function of $`\stackrel{}{p}`$ and $`t`$ alone. The kinetic equation above, in the absence of collisions, has the free streaming solution for $`f`$ , namely,

$$f(t,p)=\mathrm{exp}\left(\beta _0\sqrt{p_{}^2+p_z^2t^2/\gamma _0^2}\right),$$
(11)
where $`\beta _0`$ and $`\gamma _0`$ are constants.

We shall now consider the collision term $`C[f]`$ on the right hand side of the above equation. If we assume that the changes in momentum $`q`$ in a collision are small, namely, $`q/Q_s1`$, we can treat the collision integral as simple diffusion in momentum space. Following the discussion in Lifshitz and Pitaevskii , one can write

$`C[f]={\displaystyle \frac{s_\alpha }{p_\alpha }},`$ (12)
where the flux

$$s_\alpha =\left(\frac{\alpha N_c}{\pi }\right)^2d^3\stackrel{}{p}\left(f\frac{f^{}}{\stackrel{}{p^{}}_\beta }f^{}\frac{f}{\stackrel{}{p}_\beta }\right)\left(\delta _{\alpha \beta }(1\stackrel{}{v}\stackrel{}{v^{}})+v_\alpha v_\beta ^{}+v_\beta v_\alpha ^{}\right).$$
(13)
Note that $`\stackrel{}{v}=\stackrel{}{p}/|\stackrel{}{p}|`$.

We have used in Eq. 13 the cross–section for gluon–gluon elastic scattering, which in the limit of small momentum transfer squared $`q^2=\widehat{t}`$, is

$$\frac{d\sigma }{d\widehat{t}}=\left(\frac{\alpha _SN_c}{\pi }\right)^2\frac{(2\pi )^3}{2(N_c^21)}\frac{1}{\widehat{t}^2}.$$
(14)
The collision integral contains a logarithmic collinear divergence arising from small angle scattering. It is represented here by $``$, defined as

$$=_{q_{\mathrm{min}}}^{q_{\mathrm{max}}}\frac{dq}{q}.$$
(15)
The minimal momentum transfer is, $`q_{\mathrm{min}}=Q_s\theta _{\mathrm{min}}`$, where $`\theta _{\mathrm{min}}`$ is the minimum scattering angle. Also, $`q_{\mathrm{max}}`$ is the maximum momentum transfer. How one may determine $`q_{\mathrm{min}}`$ and $`q_{\mathrm{max}}`$ is discussed in the next section.

Having started from Eq. 6, combining Eq. 10 and Eq. 13, we arrive at the considerably simpler expression, the Landau equation

$$\frac{f}{t}\frac{p_z}{t}\frac{f}{p_z}=\lambda n_\stackrel{}{p}^2f+2\lambda n_1_\stackrel{}{p}(\stackrel{}{v}f),$$
(16)
with the definitions

$$n(t)=g_G\frac{d^3p}{(2\pi )^3}f(t,\stackrel{}{p}),n_1(t)=g_G\frac{d^3p}{(2\pi )^3}\frac{f(t,\stackrel{}{p})}{|\stackrel{}{p}|},\lambda =2\pi \alpha _S^2\frac{N_c^2}{N_c^21}.$$
(17)
Here, $`g_G=2(N_c^21)`$ and $`_\stackrel{}{p}`$ refers to differentiation with respect to $`\stackrel{}{p}`$.

Multiplying both sides of Eq. 16 by $`d^3p`$ and integrating, one finds

$`{\displaystyle \frac{d}{dt}}(tn)=0,`$ (18)
namely, the number density of gluons has the behavior $`n1/t`$. Taking the second moment of Eq. 16, one finds the following exact relation in the central region

$`{\displaystyle \frac{ϵ}{t}}+{\displaystyle \frac{\left(ϵ+P_L\right)}{t}}=0,`$ (19)
where the local energy density is given by

$`ϵ(t)=g_G{\displaystyle \frac{d^3p}{(2\pi )^3}pf(p,t)},`$ (20)
and the longitudinal pressure is

$`P_L(t)=g_G{\displaystyle \frac{d^3p}{(2\pi )^3}\frac{p_z^2}{p}f(p,t)}.`$ (21)
In the hydrodynamic regime where $`P_L=ϵ/3`$, the ideal gas equation of state gives $`ϵt^{4/3}`$. We will return to a discussion of these quantities later on in the paper.

Following Mueller, it is convenient to define the quantities $`\eta =tn`$ and $`\eta _1=tn_1`$. The Landau equation can thus be re–written as

$$t\frac{f(t,\stackrel{}{p})}{t}p_z\frac{f(t,\stackrel{}{p})}{p_z}=\lambda \eta _\stackrel{}{p}^2f(t,\stackrel{}{p})+2\lambda \eta _1(t)_\stackrel{}{p}(\stackrel{}{v}f(t,\stackrel{}{p})).$$
(22)
If $`\eta _1`$ is taken to be independent of time, this equation is the Fokker–Planck equation: the coefficient of the first term would be the diffusion coefficient, while that of the second would be the coefficient of friction.

## 4 Screening of collinear divergences in the Landau transport equation

In the one gluon exchange approximation, as in Coulomb scattering, the transport integral in Eq. 13 has a logarithmic collinear divergence–represented by $``$. In Ref. , the condition applied at early times was that in the maximal distance $`1/q_{\mathrm{min}}`$ corresponding to the minimal momentum transfer $`q_{\mathrm{min}}`$, one has at most one scattering  This condition satisfies the requirement that the soft field (due to all the other hard gluons) that a particular hard gluon scatters off be at most of size $`1/g`$. This is so because the freed gluons come from the saturated component of the nuclear wavefunction of size $`1/g`$–see the discussion in section 2.. From this condition, by using the initial gluon distribution, Mueller finds at the earliest times that

$`M_m^2=\left({\displaystyle \frac{c\alpha _SN_c}{\pi }}\right)^{2/3}{\displaystyle \frac{Q_s^2}{(Q_st)^{2/3}}}.`$ (23)
A more general expression, for any arbitrary number density , is easily derived to be

$$M_m^2=\left(\left(\frac{\alpha _SN_c}{\pi }\right)^2\frac{(2\pi )^3}{2(N_c^21)}n(t)\right)^{2/3}.$$
(24)
Alternately, one can use the more conventional prescription of regulating the logarithmic collinear divergence in $``$ through the exchange of dressed gluons characterized by a Debye screening mass. A kinetic theory expression for the Debye screening mass can be derived from the screening of a time-like gluon propagator in a medium of gluon excitations . One obtains

$$m_D^2=\frac{\alpha _SN_c}{\pi ^2}\underset{|\stackrel{}{q}|0}{lim}d^3p\frac{|\stackrel{}{p}|}{\stackrel{}{q}\stackrel{}{p}}\stackrel{}{q}\stackrel{}{}_\stackrel{}{p}f(t,\stackrel{}{p}).$$
(25)
One can check that for an equilibrium Bose distribution, one recovers the standard result: $`m_D^2=4\pi \alpha _SN_cT^2/3`$. For our initial conditions, it is the screening of the longitudinal gluons in the transverse direction that is relevant. Performing an integration by parts, this transverse Debye mass is given by the expression

$`m_T^2={\displaystyle \frac{\alpha _SN_c}{\pi ^2}}{\displaystyle \frac{d^3p}{|\stackrel{}{p}|}f(p)}.`$ (26)
At the very early stage of the $`22`$ scattering, the expression in Eq. 26 is, in principle, less reliable since the system is completely out of equilibrium. However, at later times, as the system approaches equilibrium, $`m_D`$ should be more reliable. One could attempt to parameterize the infrared cut-off of small angle scattering such that it interpolates between the two limits. However, at present there is no theoretical justification of any particular form. Since collisions become more frequent as one approaches equilibrium, it is likely that the Debye mass is a more reliable cut-off to use. At any rate, it is important that one at least have the right limit as the system approaches equilibrium.

In Refs. , the value of $`q_{\mathrm{max}}`$ is held fixed at early times. This is because if one re–writes $``$ in Eq. 15 as

$`=\mathrm{ln}\left({\displaystyle \frac{p_t}{m_T}}\right),`$ (27)
it is argued, on the basis of analytical approximations to Eq. 31 below, that $`\frac{p_t}{Q_s}`$ changes appreciably only over a long time scale ($`te^{\frac{1}{\sqrt{\alpha _S}}}`$). Since the Debye mass has a stronger time dependence, $``$ grows with time. On account of the approximations made in arriving at this conclusion, it is useful to check whether this behavior of $``$ is obtained by computing $`p_t`$ and $`m_T`$ self–consistently at each step in our numerical simulations \**\**\**We thank A. Dumitru for discussions on this point.. We indeed do find, as shown in in Fig. 1, that $``$ grows rapidly initially as a function of time but levels out to a constant as the system approaches equilibrium .

In Fig. 1, we plot the infrared masses in Eqs. 24 and 26, as well as the corresponding contribution $``$ to the collision integral, as a function of time (in units of $`1/Q_s`$). One sees that while the relative difference at early times is large, it is much less so at later times. We have also solved the transport equation using both cut-offs. While there is some quantitative difference between the two, they do not affect our qualitative conclusions. For the rest of this work, we will therefore use the Debye mass in Eq. 26.

Debye screening occurs due to the exchange of longitudinal gluons. One may worry whether the screening of transverse gluons will introduce an additional scale since the static magnetic mass is parametrically larger than the Debye mass. However, the conventional wisdom is that, for transport cross–sections, magnetic screening is also cut off in the infrared by the Debye mass . The same, for instance, is not true of color transport–the color conductivity is regulated by a magnetic mass of order $`g^2T`$ . The dynamical screening of infrared divergences in $`2n`$ ($`n2`$) transport cross-sections merits further study.

## 5 Results from numerical solution of the Fokker-Planck Equation

We will first consider, as a “warmup exercise”, the solution of Eq. 22 for the case where $`\eta _1`$ is held constant. The Landau equation then reduces to the Fokker–Planck equation. To check the consistency of our numerical simulation, we will compare, in detail, the numerical solution of the Fokker–Planck equation to the early time analytical solution found by Mueller in Refs.  and . We will discuss solutions to the non–linear Landau equations in the next section.

Following Mueller, we first define $`\stackrel{~}{f}=tf`$, $`=\frac{1}{3}\mathrm{log}(t/t_0)\frac{1}{3}\xi `$, where $`t_0`$ is a constant defined in Eq. 5. (As discussed in the previous section, this assumption is justified at early times.) Scaling $`\lambda /3\lambda `$, we can rewrite (22) as

$$\left(\frac{}{\xi }\frac{}{p_z}p_z\right)\stackrel{~}{f}=\lambda \eta \xi _\stackrel{}{p}^2\stackrel{~}{f}+2\lambda \eta _1_\stackrel{}{p}(\stackrel{}{v}\stackrel{~}{f}).$$
(28)
Eq. (28) is a Fokker–Planck equation if $`\eta _1`$ is held fixed. If one integrates Eq. 28 with respect to $`d^3p`$, $`d^3pp_{}^2`$ and $`d^3pp_z^2`$ we obtain Eqs. (17), (18) and (19) in Ref. ,

$$\frac{d}{d\xi }\eta =0,$$
(29)
$$\frac{d}{d\xi }p_z^2+2p_z^2=2\lambda \eta \xi \frac{\lambda \eta _1\xi }{\eta }\stackrel{~}{P}_L(\xi ),$$
(30)
and

$$\frac{d}{d\xi }p_{}^2=4\lambda \eta \xi \left(1\frac{\eta _1\eta _{+1}}{\eta ^2}\right)+\frac{\lambda \eta _1\xi }{\eta }\stackrel{~}{P}_L(\xi ),$$
(31)
where $`\stackrel{~}{P}_L=tP_L`$ is defined by eq. (21). These equations are exact. We can therefore perform the following consistency check on our solution. With the initial conditions in Eq. 3 and the boundary condition that $`f(p_x,p_y,p_z,t)=0`$ if any $`p_x`$, $`p_y`$ or $`p_z\mathrm{}`$, we numerically determine $`f`$ for all time, directly compute $`p_{}^2`$ and $`p_y^2`$, and check to see if Eqs. (29) to (31) are satisfied. Mueller supposes that at early times the first term on the l.h.s and the last term on the r.h.s of Eq. 30 are small, thereby yielding the solution

$$p_z^2=\lambda \eta \xi .$$
(32)
It is demonstrated in Fig. 2 that this approximation, in practice, is not a very good one. At very early times, if $`\xi 0`$ we see that the first term on the l.h.s of Eq. 30 is of the same order as the remaining terms. We cannot therefore ignore it. We can ignore the term proportional to $`\stackrel{~}{P}_L(t)`$ in Eq. (30) for early times. Then Eq. (30) becomes:

$$\frac{d}{d\xi }p_z^2+2p_z^2=2\lambda \eta \xi .$$

Its solution is

$$p_z^2=\lambda \eta \xi +\frac{\lambda \eta }{2}(\mathrm{exp}(2\xi )1).$$
(33)
This solution (see Fig. 4) agrees much better with the numerical results than Eq. (32) does. We define $`C_I(\xi )`$ as

$$2\lambda \eta \xi \frac{\lambda \eta _1\xi }{\eta }\stackrel{~}{P}_L(\xi )\frac{d}{d\xi }p_z^22p_z^2=C_I(\xi ),$$
(34)
where $`C_I(\xi )`$ should be small if Eq. (30) is satisfied. Fig. 2 verifies indeed that this is so.

Next, look at Eq. (31). If we assume that both $`\eta _1`$ and $`\eta _{+1}`$ are constant with respect to time and determined by the initial condition Eq. (3),

$$\eta =c\frac{N_c^21}{4\pi ^2\alpha _SN_c}Q_s^2,\eta _1=c\frac{N_c^21}{2\pi ^2\alpha _SN_c}Q_s,\eta _{+1}=c\frac{N_c^21}{6\pi ^2\alpha _SN_c}Q_s^3.$$
(35)
Using Eq. (35) we re–write Eq. (31) as:

$$\frac{d}{d\xi }p_{}^2=\frac{4}{3}\lambda \eta \xi +\frac{\lambda \eta _1\xi }{\eta }\stackrel{~}{P}_L(\xi ).$$
(36)
Again the term containing $`\stackrel{~}{P}_L`$ is small. We therefore solve the equation above, with the appropriate boundary conditions, to find

$$p_{}^2=\frac{Q_s^2}{2}\frac{2}{3}\lambda \eta \xi ^2.$$
(37)
As previously, define $`C_{II}(\xi )`$ such that

$$\frac{d}{d\xi }p_{}^2+\frac{4}{3}\lambda \eta \xi \frac{\lambda \eta _1\xi }{\eta }\stackrel{~}{P}_L(\xi )=C_{II}(\xi ).$$
(38)
Fig. 3 verifies that Eq. (38) is satisfied for early times.

Finally, we check to see how well Eq. (29) is satisfied numerically. Fig. 4 plots the numerical result for $`\eta `$ vs. $`t`$ and demonstrates that $`\eta `$ is constant with respect to $`t`$. As discussed earlier, Fig. 4 also compares the numerical result for $`p_z^2`$ and $`p_{}^2`$ with the analytical solutions of Eqs. (36) and (30) given by Eqs. (33) and (37).

## 6 Results from numerical solution of the Landau transport equation

In the previous section, we numerically solved the Boltzmann equation using Mueller’s approximation for the minimum scattering angle, and with $`\eta _1`$ fixed for all time (Eq. (28)) ††††††Recall that $`\eta _1`$ is an integral representing the “-1” moment of $`f`$.. Since the previous analysis is only expected to be valid for early times, one must solve the general expression (22) with $`\eta _1`$ determined self consistently.

In this section, we will discuss results from numerical solutions of the Landau equation in Eq. 22 employing the initial single particle distribution given by Eq. 3 and the dynamical screening mass given by Eq. 26. The equation is a second order partial integro–differential equation. It can be solved by combining the Alternating Direction Implicit (ADI) method and the Crank–Nicholson differencing scheme . A detailed discussion of the numerical procedure can be found in Appendix A.

We will first begin by briefly describing the evolution of the single particle distribution and its moments. Next, we will discuss various measures of kinetic equilibration and hydrodynamic flow for a system undergoing boost–invariant one dimensional expansion. Finally, we will discuss our numerical results for these quantities for different values of $`Q_s`$ and $`\alpha _S`$ ‡‡‡‡‡‡We will ignore running coupling effects in this analysis. It is argued in Ref.  that running coupling effects in the evolution are O($`\sqrt{\alpha _S}`$) and therefore suppressed..

### 6.1 Single particle distributions

The initial distribution for the numerical solution of Eq. (22) is given by (3), at the initial time, determined from Ref. , given by Eq. (4). The time evolution of the single particle distribution $`f(t,p_z,p_{})`$ is shown in Fig. 5 for a particular set of initial parameters. It begins its time evolution as a delta function in $`p_z`$ (represented practically by a narrow Gaussian distribution) and as a step function in $`p_{}`$ — see Eq. (3). As time proceeds, we see the $`p_{}`$ modes get scattered off the transverse plane in the longitudinal direction. The larger $`p_{}`$ modes decrease rapidly as the $`p_z`$ distribution widens out. At about $`t=11.95\mathrm{GeV}^1`$, the $`p_z`$ distribution is at its broadest extent.

The behavior of the single particle distributions are seen more clearly in Fig. 6, where we have plotted $`f`$, as a function of $`p_z`$ or $`p_x`$, for fixed $`p_x`$ and $`p_z`$ respectively. We see initially that the two distributions are widely different. At time proceeds, they converge relatively rapidly at soft momentum but much more slowly in the tails. If thermal equilibrium is defined strictly as the distributions being completely isotropic, this condition is reached only asymptotically, if at all.

Fig. 7 shows, for particular initial parameters, the averages $`p_x^2`$, $`p_y^2`$ and $`p_z^2`$ as functions of time. $`p_z^2`$ starts at zero and quickly rises before converging slowly to $`p_x^2`$ and $`p_y^2`$, which, in turn, decrease monotonically with time. The non–trivial behavior of $`p_z^2`$ is because the system is undergoing longitudinal expansion. In a box at rest, one expects that $`p_z^2`$ will show more of a monotonic behavior before leveling off. It again appears from Fig. 7 that the convergence to isotropic distributions is very slow. This is particularly so since the second moment of the distribution weights the high momentum tail unduly. The latter, as we observed in Fig. 6, takes longer to equilibrate.

Interestingly, as we will discuss in the following, the convergence of bulk thermodynamic observables to the expected equilibrium values is much more rapid. A likely explanation is that these observables are much less sensitive to the high momentum tail of the distribution.

### 6.2 Kinetic equilibrium and hydrodynamics

From Figs. (6) and (7), it appears that the distribution, as a whole, become isotropic only asymptotically. This statement is particularly true of the tails, in $`p_z`$ and $`p_t`$, of the distributions. The distributions agree more closely at softer momenta. Nevertheless, as we shall discuss below, when we look at thermodynamic signatures of equilibrium, they converge relatively rapidly to the expected behavior in a fluid undergoing boost–invariant 1–dimensional expansion.

Let us first discuss what this expected behavior is. Since we only consider $`22`$ processes in this work, the total number of gluons is fixed. The equilibrium solution of the classical Boltzmann equation, Eq. (22), thus has the form

$$f(t,p)=\mathrm{exp}\left[\beta (t)(\mu (t)p)\right].$$
(39)
Here $`\mu (t)`$ is the chemical potential, and $`\beta (t)=1/T(t)`$, where $`T`$ is is the temperature. Substituting this equilibrium distribution in Eq. (17), one obtains

$$n(t)=2\frac{(N_c^21)}{\pi ^2}T^3(t)\mathrm{exp}\left(\frac{\mu }{T}\right),n_1(t)=\frac{(N_c^21)}{\pi ^2}T^2(t)\mathrm{exp}\left(\frac{\mu }{T}\right).$$
(40)
The average energy density is defined to be

$$ϵ(t)=g_G\frac{d^3p}{(2\pi )^3}|\stackrel{}{p}|f(t,p).$$
(41)
From Eq. (41) and Eq. (40), the energy per particle at equilibrium is,

$$E(t)=ϵ(t)/n(t)=3T(t).$$
(42)
Now, the entropy density of a classical Boltzmann gas is defined as

$$s(t)=g_G\frac{d^3p}{(2\pi )^3}f\mathrm{log}f.$$
(43)
In equilibrium, the entropy per particle is simply,

$$S(t)=\frac{s(t)}{n(t)}=3\frac{\mu }{T}.$$
(44)
We noted previously, since the number of gluons is conserved, that $`tn(t)=constant`$. Also note that, since the entropy per particle is constant in equilibrium, we find from the above equation that $`\mu /T=constant`$. From these constraints, the system in equilibrium must satisfy

$`T^3t=\mathrm{constant},ϵt^{4/3}=\mathrm{constant}.`$ (45)
Finally, recall we had defined the longitudinal pressure in the central slice in Eq. 21 as

$`P_L(t)=g_G{\displaystyle \frac{d^3p}{(2\pi )^3}\frac{p_z^2}{p}f(p,t)}.`$
One can similarly define the transverse pressure $`P_T`$ to be

$`P_T(t)=g_G{\displaystyle \frac{d^3p}{(2\pi )^3}\frac{p_t^2}{2p}f(p,t)}.`$ (46)
From Eq. 19, and the above relations, the condition for ideal hydrodynamics is

$`P_T=P_L={\displaystyle \frac{1}{3}}ϵ.`$ (47)
The approach to equilibrium, in the sense of “saturating” the above thermodynamic (and hydrodynamic) identities, has been studied, in the relaxation time approximation, by several authors . This approximation is also employed in studies with mini–jet initial conditions . The collision kernel in Eq. 10 can be written formally as  as

$`C[f(p_z,p_t,t)]={\displaystyle \frac{\left(ff_{\mathrm{equil}}\right)}{\theta }},`$ (48)
where $`\theta `$, the collision or relaxation time, is in general a function of time and momentum. The relaxation time approximation is one where the momentum dependence of $`\theta `$ is neglected. Baym studied equilibration in this approximation taking $`\theta `$ to be a constant. He showed that the thermodynamic relations in Eqs. 45 and 47 were satisfied only asymptotically in time 111This is indeed what we would conclude from our study of the single particle distributions in the previous sub–section.. Nevertheless, the convergence of the system to the asymptotic value is rapid. The system, while not quite in local thermodynamic equilibrium, is sufficiently close to it that equilibrium is a good working assumption. Subsequently, Gavin studied equilibration in the relaxation time approximation assuming $`\theta `$ to have the time dependence ; $`\theta =\alpha t`$, where $`\alpha `$ is a constant. Depending on the value of $`\alpha `$, the system approaches the hydrodynamic limit quickly (smaller $`\alpha `$’s) or free streaming (larger $`\alpha `$’s). Heiselberg and Wang , studied the general case, $`\theta =t^p`$. They conclude that for $`p<1`$ thermalization is attained, while the system free streams for $`p>1`$. The case $`p=1`$ studied by Gavin is the marginal one, interpolating between the two regimes for different $`\alpha `$’s. In a subsequent paper , Heiselberg and Wang study the dependence of $`\theta `$ in finite temperature QCD, and tentatively conclude that $`p1/3<1`$.

### 6.3 Convergence to the hydrodynamic limit

In this work, as discussed in section 3, we go beyond the relaxation time approximation in treating the collision kernel. It is not necessary therefore, once a screening mechanism is postulated, to specify the time dependence of the collision time. In the following, we will discuss our results for the thermodynamic relations stated in the previous section.

In Fig. 8 we show the entropy per particle, $`S`$, plotted as a function of time. As shown for typical values of $`\alpha _S`$, $`c`$ and $`Q_s`$, it initially increases rapidly and flattens out, monotonically approaching its asymptotic value (thereby implying that $`\mu /T`$ goes to a constant as well). The expression in Eq. 43, used in the computation of Fig. 8, is the correct one, except at very early times when the well–known full expression for the entropy density of a Bose gas should be used. Thus, where the result in Fig. 8 is going to zero is where the classical Boltzmann expression for the entropy density is breaking down. We have checked that the full expression ensures that the entropy per particle is always positive definite. The triangles in Fig. 8 correspond to the entropy per particle assuming the partons are initially localized on the two dimensional transverse plane.

The convergence of other thermodynamic quantities to the hydrodynamic limit is shown in Fig. 9. Again, we note that the convergence to their asymptotic values is much more rapid than one would expect by looking at the single particle distributions alone. Our results suggest that the collision time effectively has a time dependence $`t^p`$ with $`p<1`$.

Since the convergence to the hydrodynamic limit is only asymptotic, deciding when the system can be described in terms of thermodynamic quantities is somewhat subjective. Here, we define the equilibrium time $`t_{eq}`$ as the time it takes for $`Tt^{1/3}`$ (42) and the entropy per particle $`S`$ (44), to reach $`90\%`$ of their maximum asymptotic value (see Figs. 8 and 9). One can then also extract the “initial” temperature and chemical potential that correspond to $`t_{eq}`$ by using Eqs. 4044.

The equilibration time $`t_{eq}`$ is very sensitive to the values of $`\alpha _S`$, $`c`$ and $`Q_s`$. Table I shows the equilibration time $`t_{eq}`$, the initial temperature $`T_{eq}`$, and the chemical potential $`\mu `$ for typical values of $`\alpha _S`$, $`c`$ and $`Q_s`$. The parametric behavior of these quantities (for a fixed value of the non–perturbative constant $`c`$) as a function of $`Q_s`$, for two different values of the coupling constant $`\alpha _S`$, is also shown in Fig. 10. Larger values of $`\alpha _S`$ and $`Q_s`$ yield smaller $`t_{eq}`$.

These results can be qualitatively understood as follows. One can show (see Eq. 37) that the equilibration time is parametrically

$`t_{eq}{\displaystyle \frac{1}{Q_s}}\mathrm{exp}\left(\sqrt{{\displaystyle \frac{2\pi }{c\alpha _SN_c}}}\right).`$
In Fig. 10, we note that $`t_{eq}`$ decreases roughly as $`1/Q_s`$. Also, it is greater for smaller $`\alpha _S`$ as one would expect from this expression. Similarly, from requiring that $`tn(t)=\mathrm{constant}`$, one finds for the initial temperature that

$`T{\displaystyle \frac{Q_s}{\alpha _S\mathrm{exp}\left(\sqrt{\frac{2\pi }{c\alpha _SN_c}}\right)}}.`$
We see that this dependence on $`Q_s`$ and $`\alpha _S`$ is confirmed in Fig. 10.

At RHIC, we expect (roughly) that $`\alpha _S0.3`$, $`c1.3`$ and $`Q_s1.0`$ GeV. For the idealized initial conditions discussed here, the corresponding time and temperature at kinetic equilibrium are $`t_{eq}3.2\mathrm{fm}`$ and $`T_{eq}174\mathrm{MeV}`$. At LHC, we expect $`Q_s2`$$`3`$ GeV, therefore $`t_{eq}2.4\mathrm{fm}`$$`1.4\mathrm{fm}`$ and $`T_{eq}321\mathrm{MeV}`$$`472\mathrm{MeV}`$. Are these numbers realistic? Not so (especially at RHIC), since it is unlikely that, at these temperatures, the system is a weakly coupled gluonic gas. At the time scales and temperatures corresponding to a rapid convergence of the system to the hydrodynamic limit, other (likely non–perturbative) effects might become important. One cannot conclude definitively that this is the case because we have not considered realistic initial distributions nor have we discussed the importance of particle number changing processes. We will comment on this point in the final section below.

## 7 Summary and outlook

We have solved numerically a non–linear transport equation, Eq. (22), which describes the evolution, after a heavy ion collision, of single particle gluon distributions in the central rapidity slice. The initial conditions for the solution of this equation are the highly anisotropic, idealized, initial conditions first discussed in the context of parton transport by Mueller . The distributions are controlled by a single scale–the saturation scale $`Q_s`$ of parton distributions in the nuclear wavefunction before the collision. Only elastic gluon–gluon scatterings are treated. Equilibration, in this approach, is dominated by small angle scattering. The collinear divergences that occur are regulated dynamically by a cut-off of the Debye form. We have checked that, at early times, the linearized Fokker–Planck equation reproduces the analytical results of Mueller. These analytical approximations however cannot be carried through to times relevant for equilibration.

We find that the tails of the initially highly anisotropic initial conditions converge very slowly to the expected isotropic equilibrium distribution. This behavior is confirmed by the behavior of $`p_z^2`$ and $`p_x^2`$. Despite the slow convergence of the single particle distributions to the isotropic thermal shape, thermodynamic observables such as the entropy per particle, the energy per particle, and the transverse and longitudinal pressures converge more rapidly to the hydrodynamic behavior expected of 1–dimensional, boost invariant expansion. The more rapid convergence for these observables occurs because they are more sensitive to softer momentum modes and less so to the high momentum tails.

Using a particular criterion for equilibration (that thermodynamic observables have reached $`90\%`$ of their asymptotic value), we extracted the equilibration times and the initial temperatures and chemical potentials. We have studied how they behave for varying values of $`Q_s`$ and $`\alpha _S`$. Even though small angle scattering is very inefficient (the equilibration time is an order of magnitude greater than the formation time), it is nevertheless smaller than the hydrodynamic time scale $`t_{hydro}R/c_0`$. Here $`R`$ is the radius of the nucleus and $`c_0`$ is the speed of sound in the fluid. The relatively long equilibration times correspond to relatively low temperatures, $`TQ_s`$. At these temperatures, it is unlikely that the system can be described as a weakly coupled gluon gas. It is therefore reasonable to ask whether other effects, be they perturbative or non–perturbative in nature, may significantly alter our results 222For an interesting recent take on this topic see Ref. .. We will enumerate below, in order of increasing complexity, those effects that are amenable to a weak coupling treatment.

Firstly, recall that our results were obtained for idealized initial conditions. More realistic initial distributions have recently become available , and they are qualitatively different from the idealized distributions. Secondly, since the occupation numbers of the gluons are large initially, the effects of Bose enhancements in the Boltzmann equation should be taken into account. Both of these effects are straightforward to incorporate in our approach, and results with these improvements will be reported shortly .

Thirdly, as we discussed in the introduction, number changing $`23`$ processes may be very important in establishing kinetic as well as chemical equilibrium . Even though these processes are suppressed by an additional power of $`\alpha _S`$, they are more efficient in re–distributing the momenta. Unfortunately, past treatments have been handicapped by uncertainties in how one treats infrared divergences. Work in progress suggests that this problem may be cured, and that equilibration is much more rapid . If this is indeed the case, then initial temperatures are relatively close to the initial saturation scale, and the weak coupling treatment of equilibration will have been self–consistent.

Ultimately, the goal is relate the rich variety of hadronic and electromagnetic spectra that will soon be available at RHIC (and some years later, at LHC) to properties of the initial nuclear wavefunction, that may independently be probed in deeply inelastic eA or in pA collisions. This paper is a first quantitative step in that direction.

## 8 Acknowledgments

We thank Jean-Paul Blaizot, Edmund Iancu, Keijo Kajantie, Alex Krasnitz, Dirk Rischke, Edward Shuryak, and Xin-Nian Wang for useful comments. In particular, we would like to thank Adrian Dumitru and Miklos Gyulassy on the one hand, and Al Mueller and Dam Son on the other, for illuminating discussions of their respective works in progress. Larry McLerran has provided wise comments and encouragement throughout the course of this work. This work was supported under DOE Contract No. DE–AC02–98CH10886 at BNL and DE-FG02-87ER40328 at the University of Minnesota.

## Appendix A Numerical analysis of second order partial integro-differential equations

Eq. (22) is a second order partial integro-differential equation whose numerical solution requires some care to ensure stability. Fortunately, standard finite differencing schemes can be used and numerical stability guaranteed. Eq. (22), is of the form:

$$\frac{u}{t}=A(t)\left(\frac{^2u}{x^2}+\frac{^2u}{y^2}+\frac{^2u}{z^2}\right)+B(t)\left(\frac{u}{x}+\frac{u}{y}+\frac{u}{z}\right),$$
(49)
where $`u=u(t,x,y,z)`$.

To discretize Eq. (49) we first define a multi-dimensional grid:

$$t=t_0+n\mathrm{\Delta }T,x=x_0+i\mathrm{\Delta },y=y_0+j\mathrm{\Delta },z=z_0+k\mathrm{\Delta },$$

where $`n=0,1,2,\mathrm{}N`$, $`i=0,1,2,\mathrm{}I`$, $`j=0,1,2,\mathrm{}J`$, and $`k=0,1,2,\mathrm{}K`$. Next we employ the Alternating Direction Implicit (ADI) Method using the Crank-Nicholson finite differencing scheme . The ADI method is especially useful in solving parabolic multidimensional equations on rectangular grids. For problems with three spatial dimensions, such as Eq. (49), the ADI method is implemented by splitting each time step of size $`\mathrm{\Delta }T`$ into three steps of size $`\mathrm{\Delta }T/3`$. At the current fractional time step, $`n+\frac{1}{3}`$ for example, only one of the spatial derivatives are evaluated and the others are evaluated at the previous time step, $`n`$. If we choose not to split the time step in such a way, the solution to Eq. (49), after discretizing, would require us to invert a large matrix of the same size as our grid. The purpose of splitting the time step in such a way is to reduce the problem in question into the solution of a tridiagonal matrix (as shown below). To illustrate this, Eq. (49) is discretized using the ADI method with the Crank-Nicholson finite differencing scheme. First discretize and evaluate the $`x`$ derivatives at the current time step $`n+\frac{1}{3}`$:

$`{\displaystyle \frac{u_{ijk}^{n+\frac{1}{3}}u_{ijk}^n}{\mathrm{\Delta }T/3}}={\displaystyle \frac{A_n}{\mathrm{\Delta }^2}}[u_{i+1jk}^{n+\frac{1}{3}}+u_{i1jk}^{n+\frac{1}{3}}2u_{ijk}^{n+\frac{1}{3}}+u_{ij+1k}^n+u_{ij1k}^n`$
$`2u_{ijk}^n+u_{ijk+1}^n+u_{ijk1}^n2u_{ijk}^n]`$
$`+\mathrm{linear}\mathrm{terms}.`$ (50)
Increment the time step by $`1/3`$ and evaluate the $`y`$ derivatives at current time step $`n+\frac{2}{3}`$:

$`{\displaystyle \frac{u_{ijk}^{n+\frac{2}{3}}u_{ijk}^{n+\frac{1}{3}}}{\mathrm{\Delta }T/3}}={\displaystyle \frac{A_n}{\mathrm{\Delta }^2}}[u_{i+1jk}^{n+\frac{1}{3}}+u_{i1jk}^{n+\frac{1}{3}}2u_{ijk}^{n+\frac{1}{3}}+u_{ij+1k}^{n+\frac{2}{3}}+u_{ij1k}^{n+\frac{2}{3}}`$
$`2u_{ijk}^{n+\frac{2}{3}}+u_{ijk+1}^{n+\frac{1}{3}}+u_{ijk1}^{n+\frac{1}{3}}2u_{ijk}^{n+\frac{1}{3}}]`$
$`+\mathrm{linear}\mathrm{terms}.`$ (51)
Increment the time step again by $`1/3`$ and evaluate the $`z`$ derivatives at the current time step $`n+1`$:

$`{\displaystyle \frac{u_{ijk}^{n+1}u_{ijk}^{n+\frac{2}{3}}}{\mathrm{\Delta }T/3}}={\displaystyle \frac{A_n}{\mathrm{\Delta }^2}}[u_{i+1jk}^{n+\frac{2}{3}}+u_{i1jk}^{n+\frac{2}{3}}2u_{ijk}^{n+\frac{2}{3}}+u_{ij+1k}^{n+\frac{2}{3}}+u_{ij1k}^{n+\frac{2}{3}}`$
$`2u_{ijk}^{n+\frac{2}{3}}+u_{ijk+1}^{n+1}+u_{ijk1}^{n+1}2u_{ijk}^{n+1}]`$
$`+\mathrm{linear}\mathrm{terms}.`$ (52)
In Eq.’s (50) to (52), “linear terms” refers to the terms in (49) which contain first order derivatives. Since the Crank-Nicholson finite differencing scheme is stable only if the terms with the highest order of derivatives are implicit, the lower order terms are allowed to be either explicit or implicit. We assume here for brevity in notation that the linear terms are explicit and evaluated at the previous time step. Therefore, we do not write the linear terms out.

The solution to Eq.’s (50) to (52) at every time step, is simply the solution of a tridiagonal matrix at every $`n`$. For example, Eq. (50) can be arranged as follows:

$$\alpha ^nu_{i1jk}^{n+\frac{1}{3}}+\beta ^nu_{ijk}^{n+\frac{1}{3}}+\gamma ^nu_{i+1jk}^{n+\frac{1}{3}}=\delta _{ijk}^n,$$
(53)
where

$$\alpha ^n=\gamma ^n=\frac{A_n\mathrm{\Delta }T}{3\mathrm{\Delta }^2},\beta ^n=2\frac{A_n\mathrm{\Delta }T}{3\mathrm{\Delta }^2},$$

and

$$\delta _{ijk}^n=\frac{A_n\mathrm{\Delta }T}{3\mathrm{\Delta }^2}\left[u_{ij+1k}^n+u_{ij1k}^n2u_{ijk}^n+u_{ijk+1}^n+u_{ijk1}^n2u_{ijk}^n\right]+\mathrm{linear}\mathrm{terms}.$$

Eq. (53) is simply an equation of the form:

$$𝐌𝐮=𝐝,$$
(54)
where $`𝐌`$ is a tridiagonal matrix. Since we know $`𝐌`$ and $`𝐝`$ at every previous time step $`n`$, we can in principal solve for $`𝐮`$ at every current time step $`n+\frac{1}{3}`$, given the boundary conditions $`u_{ijk}^0`$, $`u_{0jk}^n`$, $`u_{i0k}^n`$ and $`u_{ij0}^n`$. The inversion of sparse matrices (such as tridiagonal matrices) is usually numerically trivial.

In solving Eq. (49) one needs to specify the boundary conditions suitable to solving a second order differential equation. In this work we have specified $`u_{ijk}^0`$ by Eq. (3) and required that $`u`$ vanish at the boundary of space: $`u_{0jk}^n`$, $`u_{i0k}^n`$, $`u_{ij0}^n`$, $`u_{Ijk}^n`$, $`u_{iJk}^n`$ and $`u_{ijK}^n=0`$. Furthermore we require that the first derivatives vanish at $`x,y,z\pm \mathrm{}`$. Therefore, $`u_{Ijk}^nu_{I1jk}^n0`$, $`u_{1jk}^nu_{0jk}^n0`$, and so on for $`j`$ and $`k`$.


format: text

text:
# A new look at simple inhomogeneous chemical evolution

## 1. Introduction

The relation between galactic star formation history and the interstellar medium (ISM) is recorded in the distribution of elements we observe today. The principal signatures are: (A) Time-integrated metallicity distribution functions (MDFs), as exhibited by long-lived stars; (B) instantaneous MDFs at a given time, for example, the present-day dispersion in the inhomogeneous ISM; (C) abundance ratios of elements with varying origin and histories; and (D) age-metallicity relations. A unified understanding of these tracers is fundamental to solving the puzzle of galaxy assembly and evolution.

Over the 30 years of work on chemical evolution modeling, the constraints offered by (B), the instantaneous MDF, have been relatively neglected. However, Tinsley (1975) showed early on that chemical inhomogeneity significantly affects the other chemical signatures. Edmunds (1975) also argued that simplistic inhomogeneous evolution of overlapping supernova remnants, with no mixing, yields a present-day ISM that is orders of magnitude smoother than what is presently observed. This issue is only rarely discussed (e.g., Roy & Kunth 1995). But recently, the significance of dispersions in the chemical signatures has been reemphasized (e.g., Wyse 1995; van den Hoek & de Jong 1997). We are therefore overdue to reexamine inhomogeneous chemical evolution.

Picking up threads from the early work, I present here an analytic model that predicts the metallicity dispersion from first principles in the limiting case of no mixing. As a reference, we return to the standard, Simple model (e.g., Tinsley 1980; Pagel & Patchett 1975; Schmidt 1963), which assumes: (1) a closed system, (2) initially metal-free, (3) constant stellar initial mass function (IMF), and (4) chemical homogeneity at all times. Most previous analytic investigations of inhomogeneous evolution were tied to the Simple model, and assumed a given dispersion in metal production: Tinsley (1975) adopted a fixed metallicity dispersion and propagated this through the Simple model; while Searle (1977) and Malinie et al. (1993) considered an ensemble of regions that individually follow the Simple model such that they yield a given dispersion.

## 2. A fresh approach

A different, perhaps simpler, approach, returns to the concept of overlapping regions of contamination, as discussed qualitatively by, e.g, Edmunds (1975). We adopt the tenets, except for (4), of the Simple model. Consider an initially metal-free ISM, in which a first generation of star forming regions is randomly distributed, occupying a volume filling factor $`Q`$. We assume that these individual regions have a distribution of metal production $`f(z)`$, which is the probability density function for obtaining metallicity $`z`$ at any point affected by star formation. Note that it is not necessary to assume anything specific about the sizes of the individual regions, but only the relative total volumes of contaminated to primordial gas, given by $`Q`$.

We now consider $`n`$ subsequent generations of star formation. The individual regions also fall randomly in the ISM, but for each generation, occupy the same filling factor $`Q`$ as the first. The role of this last assumption will become apparent in §4 below. The probability that any given point in the ISM is occupied by $`j`$ overlapping regions is given by the binomial distribution in $`n`$ and $`Q`$:

$$P_j=\left(\genfrac{}{}{0pt}{}{n}{j}\right)Q^j(1Q)^{nj},1jn;$$
(1)
and the probability of a point retaining primordial metallicity is,

$$P_0=(1Q)^n.$$
(2)
For our purposes, we define the metallicity $`z`$ to be the mass of metals per unit volume; the coexisting mass density of H is assumed to be spatially uniform. In accordance with our assumption of no mixing, we assume that the metallicity at any given point is the sum of those contributed by the individual overlapping regions $`z_i`$:

$$z=\underset{i=1}{\overset{j}{}}z_i.$$
(3)
To obtain the instantaneous MDF after $`n`$ generations, we sum the distributions $`N_j(z)`$ for individual sets of regions occupied by $`j`$ objects, weighted by their $`P_j`$:

$$N_{\mathrm{ins}}(z)=\underset{j=1}{\overset{n}{}}P_jN_j(z).$$
(4)
To derive the total MDF of all objects ever created over $`n`$ generations, we must reduce the total gas mass after each generation $`k`$ by a factor $`D_k=1k\delta `$, to account for the conversion of gas into stars. We assume the reduction increment $`\delta `$ to be constant, so that the present-day gas fraction $`\mu _1=1n\delta `$. The total, time-integrated MDF is then:

$$N_{\mathrm{tot}}(z)=\frac{1}{n}\underset{k=1}{\overset{n}{}}D_{k1}N_{\mathrm{ins},\mathrm{k}}(z)=\frac{1}{n}\underset{j=1}{\overset{n}{}}\underset{k=j}{\overset{n}{}}D_{k1}P_{j,k}N_j(z),$$
(5)
where $`P_{j,k}`$ is $`P_j`$ for a given generation $`k`$. Note that this treatment for gas consumption assumes that $`P_{j,k}`$ remain unaffected by the reduction of gas.

We now assume that $`f(z)`$ remains the same for all generations. For large $`j`$, the Central Limit Theorem applies, which predicts that the summed MDF approximates a normal distribution with mean and variance given by $`j`$ times the mean $`a`$ and variance $`\sigma ^2`$ of $`f(z)`$. Thus, $`N_{\mathrm{ins}}(z)`$ is given simply by:

$$N_{\mathrm{ins}}(z)=\underset{j=1}{\overset{n}{}}P_j(2\pi j\sigma ^2)^{1/2}e^{\frac{(zja)^2}{2j\sigma ^2}},$$
(6)
where $`P_j`$ is given by equation 1. Likewise, to estimate $`N_{\mathrm{tot}}(z)`$ for large $`k`$, the binomial distribution approximates a normal distribution with mean and variance given by

$`a_\mathrm{b}`$ $`=`$ $`kQ`$
$`\sigma _\mathrm{b}^2`$ $`=`$ $`kQ(1Q).`$ (7)
Thus equation 5 may be written:

$$N_{\mathrm{tot}}(z)=\frac{1}{n}\underset{j=1}{\overset{n}{}}\underset{k=j}{\overset{n}{}}\frac{D_{k1}}{j^{1/2}2\pi \sigma _\mathrm{b}\sigma }e^{\frac{(ja_\mathrm{b})^2}{2\sigma _\mathrm{b}^2}+\frac{(zja)^2}{2j\sigma ^2}}.$$
(8)
For small $`n`$, equations 6 and 8 break down unless $`f(z)`$ describes a normal distribution.

## 3. The form of $`f(z)`$

Several considerations suggest that $`f(z)`$ should be given by a power-law. As a rough estimate, we consider the volume affected by a star-forming region to be that of the superbubble generated by its $`N_{}`$ core-collapse supernovae (SNe). Following Oey & Clarke (1997; hereafter OC97), the mechanical power $`L=N_{}E_{\mathrm{SN}}/t_\mathrm{e}`$, where $`t_\mathrm{e}=40`$ Myr is the life expectancy of $`L`$ for OB associations, and $`E_{\mathrm{SN}}=10^{51}`$ erg is the individual SN energy. The mechanical luminosity function of $`L`$ is assumed to be a power-law with index $`\beta `$ (OC97, eq. 1). If the superbubbles grow until they are confined by the ambient pressure, the standard, adiabatic shell evolution implies that their final radii will be related to $`L`$ as $`RL^{1/2}`$ (OC97, eq. 29). The total affected volume is therefore,

$$V_sR^3L^{3/2}.$$
(9)
In contrast, the total mass of metals produced by $`N_{}`$ SNe is given by,

$$M_z=m_y(t_\mathrm{e}/E_{\mathrm{SN}})L,$$
(10)
where $`m_y`$ is the mean yield of metals per SN.

Taking $`z=M_z/V_s`$, we find that $`zL^{1/2}`$ for the superbubble population, showing that the largest superbubbles generate the lowest metallicities, owing to dilution into larger volume. The MDF for superbubbles is,

$$N_s(z)dzN_s(R)dR\frac{dz}{dR}z^{3+2\beta }dzzdz,$$
(11)
where $`N_s(R)R^{12\beta }`$ is the size distribution for objects at their final radii (OC97, eq. 41). The final relation is obtained for $`\beta =2`$, a value typical, perhaps, for all galaxies (Oey & Clarke 1998). Here we see that there are larger numbers of metal-rich objects owing to the larger numbers of smaller superbubbles. Finally, the probability density for $`z`$ at any given point in the ISM is:

$$f(z)\frac{4}{3}\pi R^3N_s(z)z^{6+2\beta }.$$
(12)
For $`\beta =2`$, we obtain,

$$f(z)dz=Cz^2dz,z_{\mathrm{min}}<z<z_{\mathrm{max}}.$$
(13)
Since $`f(z)`$ is a probability density function, its integral must be unity, and $`C`$ gives the appropriate normalization. Thus, for a given generation of star formation, $`f(z)`$ is weighted toward low-metallicity regions because of their large sizes. The limiting values $`z_{\mathrm{min}}`$ and $`z_{\mathrm{max}}`$ are determined by the corresponding $`R_{\mathrm{max}}`$ and $`R_{\mathrm{min}}`$.

For a Salpeter (1955) IMF of slope –2.35 and SN progenitor masses between 8 and 120 $`\mathrm{M}_{}`$, we use $`m_y=10\mathrm{M}_{}`$, based roughly on models by Woosley & Weaver (1995). We take $`R_{\mathrm{max}}=1300`$ pc, a rather arbitrary value corresponding to the characteristic parameter $`R_\mathrm{e}`$ of OC97; and $`R_{\mathrm{min}}=25`$ pc, corresponding to individual SN remnants. These yield $`\mathrm{log}z_{\mathrm{min}}=4.7`$ and $`\mathrm{log}z_{\mathrm{max}}=3.0`$, or, taking O as a tracer of primary elements, minimum and maximum \[O/H\] of –3.0 and –1.3 relative to solar. The adopted yield of $`m_y=10\mathrm{M}_{}`$ may be somewhat high (e.g., Woosley & Weaver 1995), but the mean $`f(z)`$ may be somewhat modified depending on $`z_{\mathrm{min}}`$ and $`z_{\mathrm{max}}`$.

## 4. Results

We can now construct Monte Carlo models of the MDF for small $`n`$. We generate the components $`N_j(z)`$ by drawing from $`f(z)`$ (equation 13) $`j`$ times and summing the drawn $`z_i`$ (equation 3). This is repeated 5000 times to obtain distributions for $`N_j(z)`$. Equations 4 and 5 then provide the instantaneous and total MDFs. The histograms in Figures 1$`a`$ and $`c`$ show these models of $`N_{\mathrm{ins}}(z)`$ and $`N_{\mathrm{tot}}(z)`$ for $`n=2`$ and $`Q=0.72`$, and panels $`b`$ and $`d`$ show the same for $`n=4`$ generations. The hatched bar in panels $`a`$ and $`b`$ shows the fraction of primordial gas remaining, as determined by $`Q`$. The analytic curves show the corresponding results from equations 6 and 8. After 2 generations at $`Q=0.72`$, we can still see the component power-law distributions in the Monte Carlo models, and there is gross deviation from the analytic approximation. However, at only $`n=4`$, we can already see how the stochastic model is approaching the analytic version. In panel $`e`$ at $`n=24`$, the models are in close agreement.

In contrast to this rapid evolution, Figure 1$`f`$ demonstrates the effect of a small $`Q`$. The model for $`n=24,Q=0.12`$ is similar to that in panel $`d`$ for $`n=4,Q=0.72`$. Although these are not identical, we can see that it is the product $`nQ`$ that characterizes the evolutionary state, as a result of equation 2. Thus, the relative filling factor of contamination has the same importance as the number of contaminating generations. While this statement may seem intuitively obvious, it is worth emphasizing, since the implications are profound.

Figure 2 shows data for the outer Galactic halo MDF (Laird et al. 1988; dash-dot line), with \[Fe/H\] converted to \[O/H\] following Pagel (1989). The solid histogram shows an inhomogeneous Monte Carlo model for $`nQ=6.4`$. There is remarkable similarity between the observations and our rudimentary model, especially in the qualitative shape, having a metal-rich tail and large ($`1.0`$ dex) dispersion. The data also agree well, coincidentally, with the Simple model (dotted line; e.g., Tinsley 1980, equation 4.3) for present-day metallicity $`\mathrm{log}z_1/\mathrm{Z}_{}=0.2`$. Both models take $`\mu _1=0`$.

However, it is essential to note that the shapes of the two models result from entirely different processes. The high-metallicity turnover in the Simple model results purely from consumption of gas into stars. Figure 3 shows the Galactic bulge, a relatively evolved system, with data from Ibata & Gilmore (1995). A Simple model having $`z_1=\mathrm{Z}_{}`$ requires $`\mu _10.25`$ to resemble the data (dotted line). The dashed line shows a similar homogeneous model with no gas consumption, i.e., a linear relation with equal numbers of stars created at all metallicities ($`\mu _1=1`$). The difference between these two homogeneous models is thus entirely due to the depletion of gas into stars, which also applies to the Simple model for the halo (Figure 2).

In contrast, the high-metallicity tail of the inhomogeneous halo model represents the vestige of the $`f(z)`$ power-law; while the low-metallicity turnover results from the Central Limit Theorem, as the component distributions progress toward gaussians (cf. Figure 1). This offers an alternative interpretation of the MDF vs. the Simple Model. In addition, it reveals the importance of the high-metallicity tail in relatively unevolved systems, for probing the parent distribution $`f(z)`$.

An evolved inhomogeneous model like the Galactic bulge should approach the same limit as the homogeneous Simple model, since the metallicities contributed by $`f(z)`$ constitute a progressively smaller fraction of the current $`\overline{z}`$ of the ISM (Edmunds 1975). Figure 3 indeed shows that all the models essentially coincide in the low-metallicity tail. However, at the high-metallicity end, the form of the inhomogeneous model mainly reflects $`\sigma _\mathrm{b}`$ for the most recent contaminating generations. These become progressively narrower (equation 2), but do not truncate like the Simple model.

The inhomogeneous model matches the observed bulge MDF for $`\mu _1=0,z_1=\mathrm{Z}_{}`$; while the Simple model would require $`z_1`$ to be 10 times larger for $`\mu _1=0`$, with almost half the distribution at $`0<`$ \[O/H\] $`1`$. Thus, it is especially noteworthy that the inhomogeneous model agrees with both the Galactic halo and bulge MDFs by varying only the single parameter, nQ. The predicted shape, dispersion, and mean $`\overline{z}`$ are linked, thus the simultaneous agreement for these features is highly encouraging. It is also clear that the supersolar MDF is a vital discriminant between the models. Malinie et al. (1993) emphasize the importance of reproducing not only the low-metallicity tail, but also the high-metallicity drop-off of the MDF.

Figure 4 shows histograms of the empirical Galactic disk G-dwarf MDF for stars in solar neighborhood, with \[Fe/H\] converted to \[O/H\] as before. The solid curve shows the analytic inhomogeneous model for $`nQ=32,m_y=10\mathrm{M}_{}`$; the dashed line shows a model with a lower yield $`m_y=1\mathrm{M}_{}`$, which is thereby 10 times more evolved for the given mean metallicity. Finally, the dotted line shows the Simple model for $`z_1=\mathrm{Z}_{}`$. All models take $`\mu _1=0.25`$.

The Simple model deviates most strongly from the data, illustrating the “G-dwarf Problem.” Although the inhomogeneous models do not solve the Problem, their MDFs are slightly more peaked, thus more similar to the data. Note also that Rocha-Pinto & Maciel (1996) truncated their derived distribution below an effective \[O/H\] $`<0.6`$. The predicted \[O/H\] FWHM dispersions in the present-day, instantaneous MDF are 0.27 and 0.08 dex for the models with $`nQ=32`$ and 320, respectively. These values can be compared to the present-day observed dispersion in the solar neighborhood evidenced, for example, by the well-known –0.3 dex deviation of the Orion nebula from solar metallicity.

In summary, I present a rudimentary, one-zone model which considers overlapping areas of contamination with no mixing. This model thus represents a limiting case opposite to the homogeneous Simple model. An under-appreciated point is that the filling factor $`Q`$ of contamination is as important as the number of star-forming generations $`n`$. The product $`nQ`$ therefore constitutes a single parameter that describes a system’s evolutionary state. $`Q`$ and $`n`$ may be independent, and are roughly associated with the global star formation efficiency and age, respectively. Thus a system with low $`Q`$ can result in a present-day metal-poor ISM, but having old stars (e.g., I Zw 18); and a high $`Q`$ can yield an old, and simultaneously, metal-rich population (e.g., the Galactic bulge). By varying only the parameter $`nQ`$, the model can remarkably match both the Galactic halo and bulge metallicity distributions, and slightly improve the disk G-dwarf Problem. As expected, relatively unevolved systems are sensitive to the parent $`f(z)`$. For example, the high-metallicity tail of the Galactic halo MDF is consistent with a power-law $`f(z)`$. However, evolved systems are independent of $`f(z)`$ and show progressively decreasing dispersions. Further development of this inhomogeneous model, and investigation of additional empirical constraints, are currently underway.

I am grateful to Cathie Clarke, Mike Fall, and Rosie Wyse for detailed discussions and suggestions, and also to Annette Ferguson for early conversations. I also enjoyed discussions with many people during recent conferences.


text:
# 1 Introduction

## 1 Introduction

The dominant hadronic decay modes of the $`B`$ meson leading to open charm in the final state involve tree-level diagrams where the $`bc`$ transition leads to a charmed meson and an external $`W`$, which often emerges as a charged $`\pi `$, $`\rho `$, or $`a_1`$. In this study, we report a new measurement of the branching fractions for $`B^0D^{}\pi ^+`$ and $`B^0D^{}\rho ^+`$ 111Here, and throughout this document, we use the convention that a particular candidate state also implies the charge conjugate state is included.. The former is known to about 10% from previous measurements, while the latter has proven more difficult to determine due to the helicity of the $`\rho `$ .

## 2 The BABAR detector and dataset

The data used in this analysis were collected with the BABAR detector while operating in the PEP-II storage ring at the Stanford Linear Accelerator Center. For this analysis, we use a sample equivalent to $`5.22\pm 0.16`$ $`\text{ fb}^1`$ of integrated luminosity, collected while running on the $`\mathrm{{\rm Y}}(4S)`$ resonance. This corresponds to $`(5.93\pm 0.21)\times 10^6`$ $`B\overline{B}`$ pairs, assuming the $`\mathrm{{\rm Y}}(4S)`$ decays only to $`B`$ mesons.

The BABAR detector is described in detail elsewhere . Charged particle tracking is provided by a 5-layer silicon microstrip detector, capable of stand-alone pattern recognition, and a 40-layer cylindrical drift chamber. This is followed outward in radius by a specialized particle identification system, based on detection of Cherenkov light generated in quartz. A calorimeter consisting of thallium-doped cesium iodide crystals provides detection for photons, and particle identification for electrons. These devices sit inside the superconducting coil, which provides a 1.5$`\mathrm{T}`$ field with 3% uniformity over the tracking volume. The flux return for the magnet is finely segmented and instrumented with resistive-plate chambers to provide both muon identification and crude hadronic calorimetry for the detection of $`K_L^0`$ mesons.

## 3 $`𝑩`$ reconstruction method

Hadronic events for this study are selected by a set of simple requirements designed to minimize systematic errors in $`B`$ counting. We require more than three charged tracks forming a primary vertex within 0.5$`\mathrm{cm}`$ of the beam spot in both transverse directions, a sum of charged and neutral energy greater than 5 $`\mathrm{GeV}`$, and the normalized second Fox-Wolfram moment for the event, $`R_2`$, calculated from charged tracks and neutrals in the $`\mathrm{{\rm Y}}(4S)`$ frame, less than 0.5. Neutrals are clusters in the calorimeter with no associated tracks, more than 30$`\mathrm{MeV}`$ in energy, and a lateral moment of the shower distribution less than 0.8. The $`R_2`$ requirement is designed to reject the jet-like continuum events over the more uniformly distributed $`\mathrm{{\rm Y}}(4S)`$ decays. A fiducial requirement of $`0.410<\theta _{Lab}<2.540`$ for charged tracks and $`0.410<\theta _{Lab}<2.409`$ for neutrals, where $`\theta _{Lab}`$ is the polar angle to the beam line, is made in calculating these quantities and in the subsequent steps of the study.

$`B^0`$ candidates in the channels $`D^{}\pi ^+`$ and $`D^{}\rho ^+`$ are reconstructed using the mode $`D^{}\overline{D}{}_{}{}^{0}\pi _{}^{}`$, followed by $`\overline{D}{}_{}{}^{0}K^+\pi ^{}`$. The $`\rho ^+`$ is seen in its decay to $`\pi ^+\pi ^0`$. All charged tracks are required to originate close to the beam spot; the distance of closest approach in the transverse plane, is required to be less than 1.5 $`\mathrm{cm}`$ and along the beam line to be less than 3 $`\mathrm{cm}`$. The daughters of the $`\overline{D}^0`$ must have a transverse momentum, $`p_T`$, greater than 100 $`\mathrm{MeV}/c`$, and include at least 20 drift chamber hits. The soft pion track is only required to have $`p_T`$ greater than 70 $`\mathrm{MeV}/c`$, taking advantage of stand-alone track finding in the silicon detector to improve acceptance. No particle identification is used for this analysis. Neutral pions are formed from pairs of photons with energy greater than 30$`\mathrm{MeV}`$ and having a lateral shower moment less than 0.8. The invariant mass of the candidate must lie within $`\pm 25`$$`\mathrm{MeV}/c^2`$ of the nominal $`\pi ^0`$ mass.

$`\overline{D}^0`$ candidates are formed from $`K^+\pi ^{}`$ combinations, where the kaon and pion are required to have a minimum momentum of 200$`\mathrm{MeV}/c`$, and the invariant mass of the candidate must be within $`\pm 2.5\sigma `$ of the nominal $`\overline{D}^0`$ mass . All $`\overline{D}^0`$ candidates are required to have momentum greater than 1.3$`\mathrm{GeV}/c`$ in the $`\mathrm{{\rm Y}}(4S)`$ frame. A vertex constraint fit is performed, for which the $`\chi ^2`$ probability must be greater than 0.1%.

The reconstructed $`\overline{D}^0`$ is then combined with a soft pion having charge opposite to the kaon to form $`D^{}`$ candidates. A vertex constraint fit is applied to improve the angular resolution for the soft pion, using a fixed effective vertical beam spot size of 40$`\mu \mathrm{m}`$. In those cases where the fit converges, $`D^{}`$ candidates are selected by the requirement that $`\mathrm{\Delta }m=m(\overline{D}{}_{}{}^{0}\pi _{}^{})m(\overline{D}{}_{}{}^{0})`$ lies within $`\pm 2.5\sigma `$ of the nominal mass difference . The width is taken from a weighted average of the two-Gaussian distribution that is required to fit the $`\mathrm{\Delta }m`$ distribution.

$`B^0D^{}\pi ^+`$ candidates are reconstructed by combining a $`D^{}`$ and a $`\pi ^+`$ with momentum greater than 500$`\mathrm{MeV}/c`$. The decay $`B^0D^{}\pi ^+`$ involves a pseudoscalar initial-state particle decaying into a vector and pseudoscalar, so that the final-state $`D^{}`$ is polarized. Therefore, the helicity angle $`\theta _H`$ between the soft pion direction and the $`D^{}`$ direction in the $`D^{}`$ rest frame should be distributed as $`\mathrm{cos}^2\theta _H`$. In contrast, combinatorial background is uniform in $`\mathrm{cos}\theta _H`$. Therefore, $`B`$ meson candidates are selected with the additional requirement $`|\mathrm{cos}\theta _H|>0.4`$. This requirement removes 40% of the background and only 5% of signal events.

For the $`B^0D^{}\rho ^+`$ mode, $`\rho ^+`$ candidates are formed by combining a $`\pi ^0`$ meson and a charged pion with momentum greater than 200$`\mathrm{MeV}/c`$. We require the momentum of the $`\rho ^+`$ to be greater than 1$`\mathrm{GeV}/c`$ and the $`\pi ^+`$$`\pi ^0`$ mass to lie within 150$`\mathrm{MeV}/c^2`$ of the $`\rho ^+`$ mass. Finally, the opening angle, $`\theta _T`$, between the thrust angle of the $`B`$ candidate and the thrust axis of the remaining charged and neutrals in the event is required to satisfy $`|\mathrm{cos}\theta _T|<0.9`$ in order to remove continuum backgrounds.

In the case of a correctly reconstructed $`B`$ meson produced by the decay of an $`\mathrm{{\rm Y}}(4S)`$, within the experimental resolution, the measured sum of neutral and charged energies, $`E_m^{}`$, must be equal to the beam energy, $`E_b^{}`$, both evaluated in the $`\mathrm{{\rm Y}}(4S)`$ frame. We define $`\mathrm{\Delta }E`$ to be the difference between the measured $`B`$ candidate energy and beam energy in the $`\mathrm{{\rm Y}}(4S)`$ frame, $`\mathrm{\Delta }E=E_m^{}E_b^{}`$. The resolution on $`\sigma _{\mathrm{\Delta }E}`$ is predicted to be $`16.0\pm 0.6`$$`\mathrm{MeV}`$ for $`D^{}\pi ^+`$ and 28$`\mathrm{MeV}`$ for $`D^{}\rho ^+`$. The beam-energy substituted mass, $`m_{\mathrm{ES}}`$, is defined as $`m_{ES}^2=\left(E_b^{}\right)^2\left(_i𝒑_i\right)^2`$, where the $`𝒑_i`$ is the momentum of the $`i`$th daughter of the $`B`$ candidate. The predicted resolution in $`m_{\mathrm{ES}}`$ is typically about 2.5 $`\mathrm{MeV}/c^2`$ for $`D^{}\pi ^+`$ and 3.1 $`\mathrm{MeV}/c^2`$ for $`D^{}\rho ^+`$. This is about a factor of 10 better than the resolution in the reconstructed invariant mass. The resolution for $`m_{\mathrm{ES}}`$ is dominated by the beam energy spread rather than by the detector resolution. It is largely uncorrelated with the error on $`\mathrm{\Delta }E`$.

The variables $`\mathrm{\Delta }E`$ and $`m_{\mathrm{ES}}`$ are used to define a signal region and also sideband regions for background studies. For all modes, the region between 5.2 and 5.3$`\mathrm{GeV}/c^2`$ in $`m_{\mathrm{ES}}`$ and between $`\pm 300`$ $`\mathrm{MeV}`$ in $`\mathrm{\Delta }E`$ is used to study the $`B`$ candidates. The peak position, $`m_B^0`$, which should be the nominal $`B`$ mass, and the resolution $`\sigma _{m_{\mathrm{ES}}}`$ are extracted from the distribution of $`m_{\mathrm{ES}}`$ after requiring $`\mathrm{\Delta }E`$ be consistent with zero to within $`\pm 2.5\sigma `$. The resolution in $`\mathrm{\Delta }E`$ is extracted from the $`\mathrm{\Delta }E`$ distribution obtained by requiring that $`m_{\mathrm{ES}}`$ lie within $`\pm 2.5\sigma _{m_{\mathrm{ES}}}`$ of $`m_B^0`$.

The signal region in the two dimensional plane $`m_{\mathrm{ES}}`$ versus $`\mathrm{\Delta }E`$ is defined as a area $`\pm 2.5\sigma `$ wide centered at the nominal $`B`$ mass, $`m_B^0`$, and at $`\mathrm{\Delta }E=0`$. By staying below $`|\mathrm{\Delta }E|=m_\pi `$, we avoid correlated background from $`B`$ decays where a real final-state charged pion is either not included in the reconstruction or a random soft pion from the recoil $`B`$ is added to the observed state.

Only one candidate per event is allowed to appear in the $`m_{\mathrm{ES}}`$ versus $`\mathrm{\Delta }E`$ distribution. The criteria selected is to consider only the entry with the smallest absolute value for $`\mathrm{\Delta }E`$. The resulting two dimensional distribution of candidate events in $`\mathrm{\Delta }E`$ and $`m_{\mathrm{ES}}`$ is shown in Fig. 1.

## 4 Branching fraction measurements

The measurement of branching fractions requires an estimate of the combinatorial background in the $`m_{\mathrm{ES}}`$ distribution near the nominal $`B`$ mass. For the channel $`B^0D^{}\pi ^+`$, the $`m_{\mathrm{ES}}`$ distribution, shown in Fig. 2, is fitted with a background function  which parameterizes how the phase space approaches zero as the energy approaches $`E_b^{}`$:

$`f_{BG}(\xi )`$ $`=`$ $`N\xi \sqrt{1\xi ^2}\mathrm{exp}(\kappa (1\xi ^2))`$ (1)
where $`\xi =m_{\mathrm{ES}}/E_b^{}`$, and the normalization and the shape are determined by the parameters $`N`$ and $`\kappa `$. The signal is characterized by a Gaussian distribution with free mass, $`m_B`$, and width, $`\sigma _{m_{\mathrm{ES}}}`$. The projection of the signal as a function of $`\mathrm{\Delta }E`$, also shown in Fig. 2, is fitted using a linear background function plus a single Gaussian distribution with free mean and width, $`\sigma _{\mathrm{\Delta }E}`$. The region $`0.3<\mathrm{\Delta }E<0.13`$$`\mathrm{GeV}`$, containing feeddown from $`D^{}\rho `$ where a pion is left out of the reconstruction, is excluded from the fit. The observed width of the signal in $`m_{\mathrm{ES}}`$ is $`2.45\pm 0.18`$$`\mathrm{MeV}/c^2`$ while the $`\mathrm{\Delta }E`$ resolution is observed to be $`18.8\pm 0.9`$$`\mathrm{MeV}`$. Acceptances are calculated in terms of the fitted widths observed in data and Monte Carlo simulation.

For the channel, $`B^0D^{}\rho ^+`$ the observed $`m_{\mathrm{ES}}`$ distribution is shown in Fig. 3. This has been fitted with the same background function and signal function described above. The projection of the signal as a function of $`\mathrm{\Delta }E`$, also shown in Fig. 3, is fitted using a linear background function plus a single Gaussian distribution with free mean and width, $`\sigma _{\mathrm{\Delta }E}`$, with the region $`0.3<\mathrm{\Delta }E<0.13`$$`\mathrm{GeV}`$ once more excluded. The width of the signal in $`m_{\mathrm{ES}}`$ is $`3.5\pm 0.3`$$`\mathrm{MeV}/c^2`$ while the $`\mathrm{\Delta }E`$ resolution is observed to be $`39.5\pm 4.7`$$`\mathrm{MeV}`$.

A detailed Monte Carlo simulation is used to determine the acceptance for the signal events. Control samples are used to determine the uncertainties on crucial performance characteristics of the simulation. The fitted yields for signal events and estimated efficiencies for each channel are listed in Table 1.

A variety of sources contribute to the systematic errors on the final branching fraction results. The number of produced $`B`$ mesons is extracted from the ratio of multihadron-to-muon pair events on- and off-resonance, after application of the simple event selection criteria described above. After extrapolation from accepted to produced numbers of $`B\overline{B}`$ pairs, the estimated uncertainty is 3.6% .

The primary check on charged track efficiency is obtained by studying the probability for observing drift chamber versus silicon detector-only tracks in inclusive $`D^0K^{}\pi ^+\pi ^+\pi ^{}`$. This is compared with the rate for finding the third track in one-versus-three topology tau-pair decays. A final check is the observed multiplicity distribution in $`\mathrm{{\rm Y}}(4S)`$ events. A systematic error of 2.5% per track with $`p_T>1`$$`\mathrm{GeV}/c`$ on the overall efficiency scale is determined by comparing the three methods. The soft pion efficiency is determined from a study of the forward-backward asymmetry in observed $`D^{}`$ decays, and is also constrained to some extent by the same charm and tau-pair studies. The $`p_T`$ resolution in the Monte Carlo events is adjusted to reproduce the resolution seen in cosmic ray muons. The overall systematic uncertainty in the final result due to tracking efficiency is estimated to be 7.9% for the final states considered in this study. An additional systematic error contribution come from changes in the ratio of observed and predicted efficiencies as the selection requirements are varied within a reasonable range (3%).

In the case of the $`B^0D^{}\rho ^+`$ channel, there are additional systematic errors due to the $`\rho ^+`$ and its polarization. The reconstruction efficiency for the $`\pi ^0`$, as modeled in our Monte Carlo simulation, is verified to be accurate to within 5% . The width of the signal in the $`\mathrm{\Delta }E`$ projection is modestly larger than predicted by the Monte Carlo simulation. Therefore, we have varied the $`\mathrm{\Delta }E`$ requirement for the signal band over a wide range around the nominal 2.5$`\sigma `$, and estimate a further 10% uncertainty due to this requirement. Finally, we have extracted the branching fraction separately in the two hemispheres of the $`\rho `$ helicity angle, $`\theta _H(\rho )`$, greater and less than zero. A further 15% systematic error is assigned on the basis of the observed difference in branching fraction. The helicity distributions for the $`\rho `$ and the $`D^{}`$ are consistent with expectations based on previous measurements  and our Monte Carlo efficiency calculation, although we have not yet attempted to determine the polarization from the data.

An additional check of the result for $`B^0D^{}\pi ^+`$ is to relax the $`\mathrm{cos}\theta _H`$ requirement, and extract the fitted signal as a function of $`\mathrm{cos}\theta _H`$ instead. The distribution of the $`B`$ candidates within $`\pm 2.5\sigma _{m_{\mathrm{ES}}}`$ of the nominal $`B`$ mass and $`\pm 2.5\sigma _{\mathrm{\Delta }E}`$ of $`\mathrm{\Delta }E=0`$ exhibits no forward-backward asymmetry after acceptance corrections as can be seen in Figure 4. This confirms our understanding of the soft pion efficiency in the $`D^{}`$ decay. Overlayed on the 126 signal candidates, with an estimated background of 4 events, is a fit with $`\mathrm{cos}^2\theta _H`$, yielding a $`\chi ^2`$ for goodness-of-fit of 8.6 for 7 degrees of freedom.

The PDG compilation  of measured branching fractions for the $`D^+D^0\pi ^+`$, $`(68.3\pm 1.4)`$%, and $`D^0K^{}\pi ^+`$, $`(3.83\pm 0.09)`$%, are used in computing our final results. The measurement errors on these branching fractions are included in the systematic error on our final result. We also assume the $`\mathrm{{\rm Y}}(4S)`$ decays into $`B^0`$$`\overline{B}^0`$ pairs with a 50% fraction; no systematic error is assigned to this value.

Based on fitted yield of signal events, the estimated efficiency, and the number of produced $`B`$ mesons in our sample, the preliminary results for the branching fractions for $`B^0D^{}\pi ^+`$ and $`B^0D^{}\rho ^+`$ are $`(2.9\pm 0.3\pm 0.3)\times 10^3`$ and $`(11.2\pm 1.1\pm 2.5)\times 10^3`$ respectively. The branching fraction for $`B^0D^{}\rho ^+`$ includes all non-resonant and quasi-two-body contributions that lead to a $`\pi ^+`$$`\pi ^0`$ invariant mass in the $`\rho `$ band. However, the acceptance for non-resonant $`D^{}`$$`\pi ^+`$$`\pi ^0`$ decays is about 15% of $`D^{}\rho ^+`$ so that, combined with the known branching fraction for this mode, the non-resonant contribution to our result for $`B^0D^{}\rho ^+`$ is expected to be quite small. Both branching fraction results compare well with previous measurements and with the world average .

## 5 Summary

$`B^0`$ decays to $`D^{}`$$`\pi ^+`$ and $`D^{}\rho ^+`$ have been studied using the decay chain $`D^{}\overline{D}{}_{}{}^{0}\pi _{}^{}`$, followed by $`\overline{D}{}_{}{}^{0}K^+\pi ^{}`$. The preliminary branching fractions obtained for these channels, $`(2.9\pm 0.3\pm 0.3)\times 10^3`$ and $`(11.2\pm 1.1\pm 2.5)\times 10^3`$ respectively, are compatible with previous observations .

## 6 Acknowledgments

We are grateful for the contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. We acknowledge support from the Natural Sciences and Engineering Research Council (Canada), Institute of High Energy Physics (China), Commissariat à l’Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), Bundesministerium für Bildung und Forschung (Germany), Istituto Nazionale di Fisica Nucleare (Italy), The Research Council of Norway, Ministry of Science and Technology of the Russian Federation, Particle Physics and Astronomy Research Council (United Kingdom), the Department of Energy (US), and the National Science Foundation (US). In addition, individual support has been received from the Swiss National Foundation, the A. P. Sloan Foundation, the Research Corporation, and the Alexander von Humboldt Foundation. The visiting groups wish to thank SLAC for the support and kind hospitality extended to them.