text:
# Chapter 1 Binary Population Synthesis: Methods, Normalization, and Surprises

## Chapter 1 Binary Population Synthesis: Methods, Normalization, and Surprises

Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA vkalogera@cfa.harvard.edu

Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA; Nicolaus Copernicus Astronomical Center, 00-716 Warszawa, Poland kbelczynski@cfa.harvard.edu

### 1 INTRODUCTION

Binary systems with compact objects (neutron stars or black holes) have provided us with some of the most surprising stellar configurations discovered in galaxies: from high– and low–mass X-ray binaries, persistent or transient and mini–quasars to binary millisecond radio pulsars and double neutron stars. For decades a lot of effort has been devoted to understanding the evolutionary history of these systems, identifying the dominant formation processes and factors that determine their properties. This multi–faceted research effort has followed two main directions. One includes studies of specific observed systems with the goal to understand their characteristics and origin, and provide us with clues about the general population properties. Although it is possible to tailor an evolutionary model to reproduce the properties of an isolated system, this exercise provides little perspective on whether the putative initial conditions and subsequent tailoring are plausible. A different and more general approach is to model the evolution of an entire ensemble of primordial binaries under a common set of assumptions. Such population synthesis models (e.g., Dewey & Cordes 1987; Politano 1988; Lipunov & Postnov 1988; de Kool 1992; Kolb 1993; and many more) can be very useful in analyzing the statistical properties of the population under study and allows comparisons to observed samples of objects.

In the first part of the paper, we discuss in some detail two different methods used in binary population synthesis calculations and the reasons for choosing one or the other (§ 2). In the second part, we focus on some recent results we have obtained from studies of the formation of X-ray binaries and double compact objects, without attempting at any level to review this area of research (§ 3, 4). In our discussion of the results we focus on issues related to the absolute normalization of models and predicted formation rates, ways of constraining them and their sensitivity to model input parameters. In § 4, we report on a new class of double neutron star systems with implications for gravitational–wave detection and gamma-ray bursts.

### 2 Population Synthesis Methods

The basic goal in population synthesis calculations is to follow the evolution of an ensemble of primordial binaries through all possible evolutionary phases until the formation of the systems of interest. Typical examples of such phases include: wind mass loss, stable or unstable and conservative or non–conservative mass and angular momentum transfer phases, formation of compact objects with associated mass loss and supernova kicks, circularization, angular momentum loss through gravitational radiation and/or magnetic braking. The details of some of these phases is not well understood, but in general the choice of prescriptions and assumptions about their treatment in the models is guided by detailed hydrodynamic or stellar evolution studies. Most often for the modeling of stellar evolution and the physical processes involved, analytic approximations adequate for statistical estimates are employed.

In what follows we discuss in some more detail two very different methods that have been used so far in population synthesis calculations.

#### 2.1 Evolution of Distribution Function in Phase Space

A class of semi-analytical methods have been employed for population synthesis calculations in various flavors over the years (e.g., Kolb 1993; Politano 1996; Kalogera & Webbink 1998; Kalogera 1998). These methods are based on the idea of evolving a distribution function through phase space using Jacobian transformations. In the specific case of binaries, the “phase space” is the space of binary parameters, i.e., masses, orbital separations, and eccentricities. A distribution function describing the population of binaries at the beginning of a given evolutionary phase, $`F(M_1,M_2,A,e)`$, can be transformed into another function $`F^{}(M_1^{},M_2^{},A^{},e^{})`$ describing the binaries at the end of this phase. The transformation can be performed using Jacobians (involves partial derivatives of final binary parameters with respect to initial parameters) and provided that the functional relationships between final and initial quantities are known:

$$F(M_1^{},M_2^{},A^{},e^{})=F^{}(M_1,M_2,A,e)\left[J\left(\frac{M_1^{},M_2^{},A^{},\mathrm{e}^{}}{M_1,M_2,A,\mathrm{e}}\right)\right]^1.$$
(1.1)
In the vast majority of cases in binary evolution, these relationships are not just known but they are such that the partial derivatives can be calculated analytically! It is often the case that the evolution of a distribution function can be calculated analytically through a whole sequence of evolutionary phases (e.g., wind mass loss, common envelope, stable mass transfer, supernova explosions). Some of the mathematical details of such cases have been derived and discussed, for example, by Kalogera (1996) and Kalogera & Webbink (1998). In its general form this computational method is semi-analytical because at times it is required that integrations are performed over physical parameters that are not of interest. Typically such integrations can be performed only numerically.

These methods of phase–space evolution have great advantages in providing us with high–accuracy results, high resolution in distributions of binary parameters at low computational costs. The numerical implementation is done at the final stage of evolution, after the binary population has been evolved through a given evolutionary sequence analytically. The grid over which the distribution function is calculated is set up on the actual set of parameters and evolutionary stage of interest and the density of this grid can be chosen to be high enough so that numerical noise is not an issue. This is to be contrasted with the concept of setting up a “grid” on the initial parameters of the population and anticipating that this initial grid will be dense enough to describe well the final properties of the population. In addition to high numerical accuracy and low computational cost, the semi-analytical population synthesis method offers physical insight to how the final properties of the population depend on the properties of their progenitors and on the physical parameters that enter the description of the various evolutionary phases. This invaluable physical insight originates from the analytical derivation of the distribution functions in each evolutionary stage and it allows us to identify the dominant physical quantities that dictate the final results. Often it also allows us to derive analytical, strict limits of physical properties (e.g., Kalogera 1996; Kalogera & Lorimer 2000) or functional dependencies between parameters.

There is however one major disadvantage to the phase-space evolution method: it can be used only if the exact evolutionary sequence is known a priori. In other words, the calculation of the evolution (in practice of the Jacobian transformations) has to be “tailored” to a specific path and, for populations of systems that are formed through a multitude of evolutionary paths, the analysis has to be repeated for each one of them. Another consequence is that population studies performed solely with this method cannot lead to the discovery of new unexplored channels of evolution.

#### 2.2 Monte Carlo Techniques

Monte Carlo methods have a straightforward application in population synthesis calculations and has been used most often in the study of binary systems (e.g., Dewey & Cordes 1987; Lipunov & Postnov 1988; Portegies–Zwart & Verbunt 1996; Han, Podsiadlowski, & Eggleton 1995; and many others). The basic design relies on the idea of generating a large number of unevolved binary and single stars with physical properties generated according to pre-determined distributions. The evolution of each “Monte Carlo object” is followed through a wide variety of evolutionary stages until the formation of systems of interest (i.e., binary compact objects). The evolution modeling relies on assumptions and prescriptions for relevant physical processes, in a way very similar to phase–space evolution calculations.

Probably the most important disadvantage of Monte Carlo techniques in population synthesis studies is that of statistical accuracy of the results. Even with current computational resources, statistical accuracy can still be a problem for studies of rare populations with relatively low formation rates, such as binary compact objects and low–mass X–ray binaries. Reducing statistical noise in predicted rates of double neutron stars, for example, to less than even 10% requires a total number of primordial binaries in excess of $`10^6`$ and poses very serious computational–cost requirements (typically 100 hours on fast workstation). Consequently, achieving a satisfactory statistical accuracy for multi-dimensional distribution of rare populations over physical parameters and not just their rate needs careful examination of the models and top–level computational resources. With respect to gaining physical insight to the properties of the population of interest and what determines them, Monte Carlo techniques are less suitable than phase–space evolution methods. For example, a strict limit on physical characteristics cannot be rigorously identified in Monte Carlo simulations. Typically only statistical statements can be made and there is always the question whether a certain area of the parameter space is physically excluded or just underpopulated.

On the other hand, among the most important advantages of Monte Carlo techniques used in population studies are: (i) they are rather simple in their conceptual design and implementation, (ii) numerical calculations can be parallelized in a straightforward manner, (iii) they allow the study of a large number of evolutionary sequences and stellar populations simultaneously, and (iv) they can lead to the discovery of new, previously unappreciated or ignored formation channels for objects of interest (see § 4).

### 3 Population Synthesis of Double Compact Objects: Predicted Rates

In what follows we discuss preliminary results from a Monte Carlo population study of double compact objects, neutron stars (NS) or black holes (BH), forming through a multitude of evolutionary channels (Belczynski, Kalogera, & Bulik 2001). In particular we focus on the absolute normalization of the models and the predicted formation rates of close NS–NS, BH–NS, and BH–BH binaries. Interest in these formation rates arises from the importance of such close systems for gravitational–wave astronomy and possibly for gamma–ray bursts. The late inspiral of close binary compact objects are primary sources for the ground–based laser interferometers currently under construction (such as LIGO, VIRGO, GEO600), and predictions of detection rates depend sensitively on the formation rates of these objects. On the other hand, the question of the origin of gamma–ray bursts (GRB) and their possible connection to compact object mergers can to be addressed to some extent through a comparison of GRB and binary compact object formation rates.

To address the above questions, though, knowledge of theoretical formation rates are essential. Population synthesis models, however, often lack a good calibration and predicted merger rates available in the literature tend to cover a large range (many orders of magnitude; for a recent review see Kalogera 2001). Therefore, it is important to use a number of ways to better constrain the absolute normalization of synthesis models.

#### 3.1 Observational Constraints

For the calibration of the population synthesis models for double compact object formation presented here, we have used a number of independent constraints on stellar populations derived from observations: (i) rate of Type II supernovae thought to be core collapse events of hydrogen–rich massive stars, (ii) rate of Type Ib/c thought to be core collapse events of hydrogen–poor stars, (iii) star formation history for our Galaxy, (iv) empirical merger rates of NS–NS binaries obtained based on the observed NS–NS sample. Additional constraints can also be obtained from (i) observed frequency of helium stars in isolation and in binaries, (v) frequency and lifetime of compact object binaries with helium stars, and (vi) formation rates of X–ray binaries.

This list of constraints includes populations and events that are directly related to the formation of double compact objects but also some others that are not, but their formation can be followed (as side products) in population synthesis calculations. We use the recent determination of supernova rates of Cappellaro, Evans, & Turatto (1999). The derived rates for Type II and Ib/c supernovae for a galaxy like the Milky Way lead to a rate ratio (II/Ibc) of $`6\pm 4`$. Estimates of the Milky Way star formation rate lie in the range $`13`$ M yr-1 (Blitz 1997; Lacey & Fall 1985). Empirical estimates of the NS–NS merger rate when taking all uncertainties into account see to lie in the range $`10^610^4`$ yr-1 (Kalogera et al. 2000).

#### 3.2 Model Assumptions

Here, we give a brief description of our population synthesis code. More details about the treatment of various evolutionary processes will be presented in Belczynski, Kalogera, & Bulik (2001).

To describe the evolution of single or non–interacting binary stars (hydrogen– and helium–rich) from the zero age main sequence (ZAMS) to carbon–oxygen (CO) core formation, we employ the analytical formulae of Hurley, Pols, & Tout (2000), whose results are in good agreement with earlier stellar models (e.g., Schaller et al. 1992). To calculate masses of compact objects formed at core–collapse events, we have adopted a prescription based on the relation between CO core masses and final FeNi core masses (Woosley 1986). Our progenitor–remnant mass relation is in agreement with the results of Fryer & Kalogera (2001) based on hydrodynamic calculations of core collapse of massive stars.

Concerning the evolution of interacting binaries, we model the changes of mass and orbital parameters (separation and eccentricity) taking into account mass and angular momentum transfer between the stars or loss from the system during Roche–lobe overflow, tidal circularization, rejuvenation of stars due to mass accretion, wind mass loss from massive and/or evolved stars, dynamically unstable mass transfer episodes leading to common–envelope (CE) evolution and spiral–in of the stars. We extend the usual treatment of CE evolution based on energy considerations (Webbink 1984) to include cases where both stars have reached the giant branch and have convective envelopes (hydrogen or low–mass helium stars). As suggested by Brown (1995) for hydrogen–rich stars, we expect the two cores to spiral–in until a merger occurs or the combined stellar envelopes are ejected. We also account for the possibility that compact objects accrete mass during CE phases (following Brown 1995). At NS formation, we model the effects of asymmetric supernovae (SN) on binaries, i.e., mass loss and natal kicks, for both circular and eccentric orbits (e.g., Kalogera 1996; Portegies–Zwart & Verbunt 1996). We assume that kicks are isotropic with a given magnitude distribution.

In the synthesis calculations, we evolve a population of primordial binaries and single stars through a large number of evolutionary stages, until coalescing double compact objects are formed (merger times $`<10`$ Gyr). The total number of binaries (typically a few million) in each simulation is determined by the requirement that the statistical (Poisson) fractional errors ($`1/\sqrt{N}`$) of the final populations are lower than 10%.

In our standard model, the properties of primordial binaries follow certain assumed distributions: for primary masses ($`5100`$ M), $`M_1^{2.7}dM_1`$; for mass ratios ($`0<q<1`$), $`dq`$; for orbital separations (from a minimum, so both ZAMS stars fit within their Roche lobes, up to $`10^5`$ R), $`dA/A`$; for eccentricities, $`2e`$. Each of the models is also characterized by a set of assumptions, which, for our standard model, are:

* Kick velocities. We use a weighted sum of two Maxwellian distributions with $`\sigma =175`$ km s-1 (80%) and $`\sigma =700`$ km s-1 (20%) (Cordes & Chernoff 1997); Black–hole kicks are weighted based on the amount of fallback mass during their formation.
* Maximum NS mass. We adopt a conservative value of $`M_{\mathrm{m}ax}=3`$ M (e.g., Kalogera & Baym 1996). It affects the relative fractions of NS and BH and the outcome of NS hyper–critical accretion in CE phases.
* Common envelope efficiency. We assume $`\alpha _{\mathrm{C}E}\times \lambda =1.0`$, where $`\alpha `$ is the efficiency with which orbital energy is used to unbind the stellar envelope, and $`\lambda `$ is a measure of the central concentration of the giant.
* Non–conservative mass transfer. In cases of dynamically stable mass transfer between non–degenerate stars, we allow for mass and angular momentum loss from the binary (see Podsiadlowski, Joss, & Hsu 1992), assuming that half of the mass lost from the donor is also lost from the system ($`1f_\mathrm{a}=0.5`$) with specific angular momentum equal to $`\beta 2\pi `$A2/P ($`\beta =1`$).
* Star formation history. We assume that star formation has been continuous in the disk of our Galaxy for the last 10 Gyr (e.g., Gilmore 2001).

An extensive parameter study is essential in assessing the robustness of population synthesis results. Apart from our standard case, we examine the results for 25 additional models, where we vary all of the above parameters within reasonable ranges. The complete set of models and the assumptions that are different from our standard choices are:

* A: standard model
* B1–7: zero kicks, single Maxwellian with $`\sigma =50,100,200,300,`$
$`400`$ km s-1, “Paczynski” kicks with $`\sigma =600`$ km s-1
* C: no hyper–critical accretion onto NS and BH in CEs
* D: maximum NS mass: $`M_{\mathrm{m}ax}=2,1.5`$ M
* E1–6: $`\alpha _{\mathrm{C}E}\times \lambda =0.1,0.25,0.5,0.75,2,3`$
* F1–4: mass fraction accreted: f$`{}_{\mathrm{a}}{}^{}=0.1,0.25,0.75,1`$
* G1–2: specific angular mom. $`\beta =0.5,2`$
* H: primary mass: $`M_1^{2.35}`$
* I1–2: mass ratio: $`q^{2.7}`$, $`q^{3.0}`$

#### 3.3 Most Important Model Parameters and their Effects

In assessing the sensitivity of our results to the model input parameters varied as described in § 3.2, we examine the behavior of (i) the calculated ratio of Type II to Type Ib/c rates, and (ii) the formation rates of three types of coalescing double compact objects: NS–NS, BH–NS, and BH–BH.

We find that the SN rate ratio is by far most sensitive to the assumed binary fraction for the stellar population. In agreement with earlier results (de Donder & Vanbeveren 1998), the ratio increases with decreasing binary fraction, primarily because the majority of Type Ib/c supernovae occur in binary systems, where loss of the hydrogen envelope is facilitated. This behavior tentatively points to a low binary fraction for our Galaxy (about 25%), or conversely to the galaxies used to determine the SN rates having a binary fraction lower than that in the Milky Way. We should note however, that the predicted SN rate ratios are consistent with the observational estimates within the one-sigma error bar.

Variations of the mass ratio distribution appear to be the most important factor for the predicted coalescence rates. A distribution that strongly favors high mass ratios leads to an increase of the NS–NS rate by a factor of about 30 relative to the rate in the case that low mass ratios are strongly favored. The BH–NS and BH–BH rates are less sensitive and show variations by factors of about 3 and 6, respectively.

The common–envelope efficiency is well known to be an important model parameter in population synthesis. Rates of NS binaries decrease by about an order of magnitude, for a similar decrease of the product $`\alpha _{\mathrm{C}E}\lambda `$, which leads to an increased rate of merger prior to the formation of the binaries of interest. However, the effect on the coalescence rate of BH–BH binaries is non-monotonic. The reason is that these binaries tend to form with wide orbits. Decreasing the CE efficiency initially causes the rate to increase, as more tight binaries are allowed to form. With further decrease of the efficiency though, the early merger rate dominates and the formation rate of BH–BH binaries drops.

Last, the effect of kicks (imparted to NS but also to BH with smaller magnitudes) affects the double compact object formation in many different ways: supernova survival, formation of tight binaries, mergers. The formation typically peak for $`\sigma =100200`$ km s-1 and decrease with increasing average kick magnitudes. The rate of close BH–NS systems is found to vary the most in our parameter study by a factor of about 3.

### 4 New Class of Close NS–NS: Without NS Recycling

Empirical rate estimates of NS–NS mergers, are derived under an important underlying assumption that all NS–NS binaries have contained a rather long–lived radio pulsar at some point in their lifetime (see Kalogera et al. 2000). However, in our very recent Monte Carlo population synthesis calculations (Belczynski & Kalogera 2000), we discovered a new evolutionary path leading to NS–NS formation that does not involve recycling of any of the two neutron stars. The implication of these findings is that there may be a significant population of old NS–NS binaries that could not be detected as binary pulsars. We use the calculated formation rate of non–recycled NS–NS relative to that of recycled NS–NS systems to derive upward revisions for the rates estimates based on the current observed sample.

#### 4.1 Model Calculations and Results

We use our population synthesis models to investigate all possible formation channels of NS–NS binaries realized in the simulations. We find that a significant fraction of coalescing NS–NS systems are formed through a new, previously not identified evolutionary path.

In Figure 1 we describe in detail the formation of a typical NS–NS binary through this new channel. The evolution begins with two phases of Roche–lobe overflow. The first, from the primary to the secondary, involves non–conservative but dynamically stable mass transfer (stage II) and ends when the hydrogen envelope is consumed. The second, from the initial secondary to the helium core of the initial primary, involves dynamically unstable mass transfer, i.e., CE evolution (stage IV). The post–CE binary consists of two bare helium stars of relatively low masses. As they evolve through core and shell helium burning, the two stars acquire “giant–like” structures, with developed CO cores and convective envelopes (e.g., Habets 1987). Their radial expansion eventually brings them into contact and the system evolves through a double CE phase (stage VI; similar to Brown (1995) for hydrogen–rich stars). During this double CE phase, the combined helium envelopes are ejected at the expense of orbital energy. The tight, post–CE system consists of two CO cores, which eventually end their lives as Type Ic supernovae. The survival probability after the two supernovae is quite high, given the tight orbit before the explosions. The end product in this example is a close NS–NS with a merger time of $`5`$ Myr (typical merger times are found in the range $`10^410^8`$ yr).

The unique qualitative characteristic of this NS–NS formation path is that both NS have avoided recycling. The NS progenitors have lost both their hydrogen and helium envelopes prior to the two supernovae, so no accretion from winds or Roche–lobe overflow is possible after NS formation. Consequently, these systems are detectable as radio pulsars only for a time ($`10^6`$ yr) much shorter than recycled NS–NS pulsar lifetimes ($`10^810^{10}`$ yr in the observed sample). Such short lifetimes are of course consistent with the number of NS–NS binaries detected so far and the absence of any non–recycled pulsars among them.

Given the uncertain absolute normalization of population synthesis models, we focus primarily on the formation rate of non–recycled NS–NS binaries relative to that of recycled pulsars, formed through other, qualitatively different evolutionary paths. Based on this comparison, for each of our models, we derive a correction factor for empirical estimates of the Galactic NS–NS coalescence rate. Since these estimates are derived based on the observed sample, they can account only for NS–NS systems with recycled pulsars, and they must be increased to include any non–recycled systems formed. For the majority of the examined models, we find this correction factor to lie in the range $`1.53`$ although for some models it can be as high as 10 or even higher (for more details see Belczynski & Kalogera 2000).

As already mentioned, the realization of the newly identified path through a double helium CE phase depends on the final stages of a helium star evolution. It has long been known that low mass helium stars, after core helium exhaustion, expand significantly and develop a “giant-like” structure with a clearly defined core and a convective envelope (Habets 1987; Avila-Reese 1993; Woosley, Langer, & Weaver 1995; Hurley et al. 2000). We further examined in detail models of evolved helium stars (Woosley 1997, private communication) and found that helium stars below 4.0 M have deep convective envelopes and that slightly more massive helium stars ($``$ 4–4.5 M) still form convective envelopes although shallower. Evolved stars with convective envelopes, overfilling Roche lobes in binary systems, transfer mass on a dynamical time scale, and as a consequence CE evolution ensues. The development of CE phases was proposed first in the context of cataclysmic variable formation (Paczynski 1976) and is now supported by detailed hydrodynamic calculations in a variety of binary configurations (e.g. Rasio & Livio 1996; Taam & Sandquist 2000 and references therein). At present no hydro calculations exist for the case of two evolved stars. Based on our basic understanding of CE development, it seems reasonable to expect that, if two stars with convective envelopes are involved in a mass transfer episode, a double core spiral-in can occur leading to double CE ejection (Brown 1995). Based on these earlier calculations, we adopt a maximum helium–star mass for CE evolution (double or single) of 4.5 M. The formation rates of both types of NS–NS binaries are somewhat sensitive to this value, because they depend on whether helium stars evolve through CE phases (single and double, for recycled and non-recycled systems, respectively). Reducing the maximum mass to 4 M actually increases the rate correction factor, although by less than 20%, as it reduces the recycled NS–NS rate by a factor larger than the non–recycled rate.

We note that the identification of the formation path for non–recycled NS–NS binaries stems entirely from accounting for the evolution of helium stars and for the possibility of double CE phases, both of which have typically been ignored in previous calculations (with the exception of Fryer et al. 1999, although formation paths for non–recycled NS–NS were not discussed).

#### 4.2 DISCUSSION

We have identified a new possible evolutionary path leading to the formation of close NS–NS binaries, with the unique characteristic that both NS have avoided recycling by accretion. The realization of this path is related to the evolution of helium–rich stars, and particularly to the radial expansion and development of convective envelopes (on the giant branch) of low–mass ($`<4.5`$ M) helium stars. We find that a significant fraction of coalescing NS–NS form through this new channel, for a very wide range of model parameters. In some cases, the non–recycled NS–NS systems strongly dominate the total close NS–NS population.

Since both NS are non–recycled, their pulsar lifetimes are too short (by $`10^3`$), and hence their detection probability is negligible relative to recycled NS–NS. However, intermediate progenitors of these systems may provide evidence in support of the evolutionary sequence. Examples are close binaries with two low–mass helium stars or a helium star with an O,B companion. Although there are selection effects against their detection too (e.g., short lifetimes, high–mass helium stars are brighter, broad spectral lines due to winds, high luminosity contrast between binary members, etc.), it has been estimated that so far only about 1% of binary helium stars have been detected (Vrancken et al. 1991). Therefore, it is reasonable to expect that the observed sample will increase in future years as observational techniques improve. It is worth noting that theoretical calculations indicate that explosions of low–mass helium stars formed in binaries reproduce the light curves of type Ib supernovae (Shigeyama et al. 1990).

The results of our population modeling show that, if one accounts for the formation of non–recycled NS–NS binaries, the total number of coalescing NS-NS systems could be higher by factors of at least 50%, and up to 10 or even higher. Such an increase has important implications for prospects of gravitational wave detection by ground–based interferometers. Using the results of Kalogera et al. (2000) on the empirical NS–NS coalescence rate, we find that their most optimistic prediction for the LIGO I detection rate could be raised to at least 1 event per 2–3 years, and their most pessimistic LIGO II detection rate could be raised to 3–4 events per year or even higher.

Our results also have important implications for gamma–ray bursts (GRBs), if they are associated with NS–NS coalescence. We find that the typical merger times of the non–recycled NS–NS are considerably shorter than those of recycled binaries, whereas their center–of–mass velocities are higher. The balance of these two competing effects could alter the current consensus for the location of GRB progenitors relative to their host galaxies (e.g., Belczynski, Bulik, & Rudak 2000; Bloom, Kulkarni, & Djorgovski 2000; Fryer et al. 1999).

### Acknowledgments

Support is acknowledged by the Smithsonian Institution through a CfA Predoctoral Fellowship to KB and a Clay Fellowship to VK, and by a Polish Nat. Res. Comm. (KBN) grant 2P03D02219 to KB.

### References


format: text

text:
# I INTRODUCTION

## I INTRODUCTION

Perturbative QCD (PQCD) factorization theorem for exclusive processes  states that a hadronic form factor at large momentum transfer is expressed as convolution of a hard amplitude with hadron distribution amplitudes. However, the PQCD evaluation of the pion form factor suffers the soft enhancement from the end point of a momentum fraction $`x0`$ . In this region the hard amplitude is characterized by a low scale, such that perturbative expansion in terms of a large coupling constant $`\alpha _s`$ is not self-consistent. More serious end-point singularities, logarithmic and linear, have been observed in the leading-twist (twist-2) and next-to-leading-twist (twist-3) contributions to the $`B\pi `$ transition form factors , respectively.

We argue that when the end-point region is important, the double logarithms $`\alpha _s\mathrm{ln}^2x`$ from radiative corrections should be organized to all orders in order to improve perturbative expansion. These double logarithms have been found in the semileptonic decay $`B\pi l\overline{\nu }`$  and in the radiative decay $`B\gamma l\overline{\nu }`$ , and resummed. Here we shall give a systematic treatment of this type of double logarithms, which is of the collinear origin, by introducing a quark jet function into PQCD factorization theorem. The all-order factorization of the jet function from the decay $`B\gamma l\overline{\nu }`$ is proved following the procedures proposed in , which provides a solid theoretical ground for the modified formalism appropriate for the end-point region. It will be shown that the jet function is defined as a matrix element of a quark field attached by a Wilson line.

The double logarithms in the jet function are resummed into a Sudakov form factor in the Mellin space using the technique developed in . In principle, with the definition of the jet function constructed in this work, its threshold resummation can be performed up to the next-to-leading-logarithm accuracy. Since our purpose is to understand the qualitative behavior of the jet function at small $`x`$, it suffices to derive only the leading-logarithm result. Extending the above formalism to the semileptonic decay $`B\pi l\overline{\nu }`$, we obtain the identical Sudakov factor, implying its universality. By means of the inverse Mellin transformation, the Sudakov factor is found to vanish quickly as $`x0`$. The suppression is so strong that the end-point singularities in the $`B\pi `$ form factors are smeared. We conclude that in a self-consistent perturbative evaluation of heavy-to-light transition form factors, the end-point singularities do not exist.

It has been argued that the dependence on the transverse momentum $`k_T`$ carried by the light spectator in a $`B`$ meson can not be dropped from hard amplitudes for exclusive $`B`$ meson decays . As considering higher-order corrections to a hard amplitude, $`k_T`$ appears in the ratio $`k_T/k^+`$, $`k^+`$ being the longitudinal component of the light spectator momentum. Obviously, this ratio is not power-like, because $`k_T`$ and $`k^+`$ are of the same order of magnitude. Hence, a reliable formalism for exclusive $`B`$ meson decays must involve the parton transverse degrees of freedom. The effect due to the inclusion of $`k_T`$ has been explored thoroughly in , which also smears the end-point singularities in the heavy-to-light form factors . In this paper we shall concentrate only on the Sudakov effect from threshold resummation, and neglect the $`k_T`$ dependence.

In Sec. II we work out the $`O(\alpha _s)`$ factorization formula for the $`B\gamma l\overline{\nu }`$ decay amplitude at the end point. The all-order factorization is proved in Sec. III, leading to the convolution of the hard amplitude with the $`B`$ meson distribution amplitudes and with the jet function. Threshold resummation is done in Sec. IV by solving an evolution equation for the jet function. In Sec. V we extract the jet function from the semileptonic decay $`B\pi l\overline{\nu }`$. The smearing effect from the Sudakov factor on the end-point singularity is demonstrated numerically by evaluating the leading-twist contribution to the $`B\pi `$ form factor. Section VI is the conclusion. In Appendix A we comment on the Sudakov resummation performed in .

## II $`O(\alpha _s)`$ FACTORIZATION

We study the radiative decay $`B\gamma l\overline{\nu }`$. The $`B`$ meson carries the momentum $`P_1=(M_B/\sqrt{2})(1,1,\mathrm{𝟎}_T)`$ and the outgoing photon carries the momentum $`P_2=(M_B/\sqrt{2})(0,\eta ,\mathrm{𝟎}_T)`$ with the energy fraction $`\eta `$, where the light-cone variables have been adopted. Consider the kinematic region with small $`q^2=(1\eta )M_B^2`$, $`q=P_1P_2`$ being the lepton pair momentum, i.e., with large $`\eta `$, where PQCD is applicable. The mass difference between the $`B`$ meson and the $`b`$ quark, $`\overline{\mathrm{\Lambda }}=M_Bm_b`$, is treated as a small scale in the heavy quark limit. As stated in , there are two types of infrared divergences in radiative corrections, soft and collinear. Soft divergences come from the region of a loop momentum $`l`$, where all its components vanish. Collinear divergences are associated with a massless parton carrying a momentum of order $`M_B`$. The $`O(\alpha _s)`$ factorization of the soft divergences into the $`B`$ meson distribution amplitudes has been presented in . Here we derive only the $`O(\alpha _s)`$ factorization of the collinear divergences at the end point into the jet function.

According to leading-twist PQCD factorization theorem, the form factor relevent to the decay $`B\gamma l\overline{\nu }`$ is written as

$`G(q^2)={\displaystyle \underset{m=+,}{}}\varphi _m(x)H_m(x,\eta ),`$ (1)
where the symble $``$ denotes the convolution over the spectator momentum fraction $`x=k^+/P_1^+`$, $`k`$ being the spectator quark momentum. The light-cone $`B`$ meson distribution amplitudes $`\varphi _m`$ are defined by

$`\varphi _\pm (x)`$ $`=`$ $`{\displaystyle \frac{dy^{}}{2\pi }e^{ixP_1^+y^{}}0|\overline{q}(y^{})\gamma _5(\mathit{}+I)\gamma ^\pm \mathrm{exp}\left[ig_0^y^{}𝑑zn_{}A(zn_{})\right]b_v(0)|B(P_1)},`$ (2)
with the dimensionless vectors $`v=P_1/M_B`$ and $`n_{}=(0,1,\mathrm{𝟎}_T)`$, and the rescaled $`b`$ quark field $`b_v`$. The hard amplitudes $`H_m`$ are obtained by contracting the quark-level diagrams with the spin structures corresponding to $`\varphi _m`$ . Figure 1 is the lowest-order example, where the upper line represents a $`b`$ quark and $`\times `$ represents a weak decay vertex. The contributions to the hard amplitudes from Figs. 1(a) and 1(b) scale like $`1/(\overline{\mathrm{\Lambda }}M_B)`$ and $`1/M_B^2`$, respectively. Below we shall concentrate on Fig. 1(a).

The factorization formula in Eq. (1) is appropriate for the region with $`k^+O(\overline{\mathrm{\Lambda }})`$, in which the only infrared divergences are the soft ones absorbed into $`\varphi _\pm `$ in Eq. (2). Near the end point $`k^+O(\overline{\mathrm{\Lambda }}^2/M_B)`$, the internal quark in Fig. 1(a) carries a large momentum $`P_2k`$ with its invariant mass vanishing like $`(P_2k)^2=2xP_1P_2O(\overline{\mathrm{\Lambda }}^2)`$. This kinematics is similar to the threshold region of deeply inelastic scattering (DIS) with the Bjorken variable $`x_B1`$, where the scattered quark also carries a momentum with a large component and possesses a small invariant mass $`(1x_B)s`$, $`s`$ being the center-of-mass energy. In this region the scattered quark produces a jet of particles, to which the radiative corrections contain additional collinear divergences. Hence, a jet function needs to be introduced into the standard factorization theorem for DIS . Similarly, a jet function has been incorporated into the factorization of direct photon production at a large photon transverse momentum (threshold) . Here we associate a jet function with the internal quark at the end point of the momentum fraction involved in the decay $`B\gamma l\overline{\nu }`$.

An additional collinear divergence from the loop momentum parallel to $`P_2`$ appears in the higher-order correction to the weak decay vertex shown in Fig. 1(c). This divergence can be extracted by replacing the $`b`$ quark line by an eikonal line in the direction of $`n_+`$:

$`{\displaystyle \frac{\mathit{}_1\mathit{}+\mathit{}+m_b}{(P_1k+l)^2m_b^2}}\gamma ^\beta b(P_1k){\displaystyle \frac{n_+^\beta }{n_+l}}b(P_1k),`$ (3)
where $`b(P_1k)`$ represents the $`b`$ quark spinor with the momentum $`P_1k`$ and $`n_+=(1,0,\mathrm{𝟎}_T)`$ is a dimensionless vector. The above replacement makes sense, since $`\gamma ^\beta `$ is dominated by its plus component $`\gamma ^+`$, $`l^{}`$ is much larger than $`l^+`$, and $`\mathit{}\gamma ^\beta `$ diminishes. The factorization of the fermion flow is achieved by inserting the Fierz identity,

$`I_{ij}I_{lk}={\displaystyle \frac{1}{4}}I_{ik}I_{lj}+{\displaystyle \frac{1}{4}}(\gamma _5)_{ik}(\gamma _5)_{lj}+{\displaystyle \frac{1}{4}}(\gamma _\alpha )_{ik}(\gamma ^\alpha )_{lj}+{\displaystyle \frac{1}{4}}(\gamma _5\gamma _\alpha )_{ik}(\gamma ^\alpha \gamma _5)_{lj}+{\displaystyle \frac{1}{8}}(\sigma _{\alpha \beta })_{ik}(\sigma ^{\alpha \beta })_{lj},`$ (4)
in which the first and last terms contribute in the combined structure,

$`I_{ij}I_{lk}{\displaystyle \frac{1}{4}}I_{ik}(\mathit{}_+\mathit{}_{})_{lj}.`$ (5)
Assigning the identity matrix $`I`$ in Eq. (5) to the trace for the hard amplitude, we obtain Fig. 1(a). The matrix $`\mathit{}_+\mathit{}_{}/4`$ then leads to the loop integral,

$`J_{}^{(1)}`$ $`=`$ $`ig^2C_F{\displaystyle \frac{d^4l}{(2\pi )^4}\frac{1}{4}tr\left[\mathit{}_+\mathit{}_{}\gamma _\beta \frac{\mathit{}_2\mathit{}+\mathit{}}{(P_2k+l)^2}\right]\frac{n_+^\beta }{n_+ll^2}},`$ (6)
$`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}C_F\mathrm{ln}^2x+\mathrm{},`$ (7)
where $`C_F=4/3`$ is a color factor. To obtain the above double logarithm, we have regularized the collinear pole from $`l`$ parallel to $`n_+`$ by allowing $`n_+`$ to possess a small amount of minus component. The correction to the photon vertex in Fig. 1(d) contains only the single logarithm $`\alpha _s\mathrm{ln}x`$, since the phase space of the loop momentum is restricted to $`0<l^+<k^+O(\overline{\mathrm{\Lambda }}^2/M_B)`$. The self-energy correction in Fig. 1(e), generating also the single logarithm, is factorized into $`J_{}^{(1)}`$ trivially by applying the Fierz transformation in Eq. (5). The explicit expression of $`J_{}^{(1)}`$ is not essential here. The reason for labeling the $`O(\alpha _s)`$ jet functions by $``$ and $``$ will become clear in the next section. To group Figs. 1(c) and 1(e), we introduce a quark jet function with the zeroth-order expression and the $`O(\alpha _s)`$ expression,

$`J^{(0)}=tr(\mathit{}_+\mathit{}_{})/4=1,J^{(1)}=J_{}^{(1)}+J_{}^{(1)},`$ (8)
respectively.

It has been shown that the soft divergences in Figs. 1(c), 1(d) and 1(f)-1(h) are absorbed into the $`O(\alpha _s)`$ $`B`$ meson distribution amplitudes $`\varphi _\pm ^{(1)}`$ . The remaining infrared finite $`O(\alpha _s)`$ contributions, including the single logarithms except that from Fig. 1(e), are assigned into the hard amplitudes $`H_\pm ^{(1)}`$. Therefore, the modified factorization in the end-point region is written, up to $`O(\alpha _s)`$, as

$`G^{(0)}+G^{(1)}={\displaystyle \underset{m=+,}{}}[1+\varphi _m^{(1)}][H_m^{(0)}+H_m^{(1)}][1+J^{(1)}]+O(\alpha _s^2),`$ (9)
with $`G^{(0)}`$ and $`G^{(1)}`$ displayed in Fig. 1. Note that $`J`$ does not depend on the subscript $`m`$, a fact obvious from the derivation of Eq. (7).

## III ALL-ORDER FACTORIZATION

In this section we prove the factorization theorem for the radiative decay $`B\gamma l\overline{\nu }`$ at the end point to all orders, and construct the definition of the jet function $`J(x)`$ via

$`J(x)\overline{u}(P_2k)u(P_2k)|\overline{q}(0){\displaystyle \frac{1}{4}}\mathit{}_+\mathit{}_{}\mathrm{exp}\left[ig{\displaystyle _{\mathrm{}}^0}𝑑zn_+A(zn_+)\right]|0.`$ (10)
The spinor $`u(P_2k)`$ is associated with the internal quark, through which the momeutm $`P_2k`$ flows. The idea of the proof is based on induction  with the help of the Ward identity,

$`l_\mu G^\mu (l,k_1,k_2,\mathrm{},k_p)=0,`$ (11)
where $`G^\mu `$ is a physical amplitude with an external gluon carrying the momentum $`l`$ and with $`p`$ external quarks carrying the momenta $`k_1`$, $`k_2`$, $`\mathrm{}`$, $`k_p`$. All these external particles are on mass shell.

The $`O(\alpha _s)`$ factorization of the collinear divergences has been worked out in Sec. II. Assume that the factorization theorem holds up to $`O(\alpha _s^N)`$:

$`G={\displaystyle \underset{m=+,}{}}\varphi _mH_mJ,`$ (12)
with the definitions,

$`G={\displaystyle \underset{i=0}{\overset{N}{}}}G^{(i)},\varphi _m={\displaystyle \underset{i=0}{\overset{N}{}}}\varphi _m^{(i)},H_m={\displaystyle \underset{i=0}{\overset{N}{}}}H_m^{(i)},J={\displaystyle \underset{i=0}{\overset{N}{}}}J^{(i)},`$ (13)
$`G^{(i)}`$ denotes the full diagrams of $`O(\alpha _s^i)`$. The $`B`$ meson distribution amplitudes $`\varphi _m^{(i)}`$ and the jet function $`J^{(i)}`$ are defined by the $`O(\alpha _s^i)`$ terms in the perturbative expansions of Eq. (2) and of Eq. (10), respectively. The $`O(\alpha _s)`$ hard amplitudes $`H_m^{(i)}`$ do not contain the intial-state soft divergences and the end-point double logarithms, which have been collected into the $`B`$ meson distribution amplitudes and into the jet function, respectively. We then have the relations,

$`G^{(k)}={\displaystyle \underset{m=+,}{}}{\displaystyle \underset{i=0}{\overset{k}{}}}{\displaystyle \underset{j=0}{\overset{ki}{}}}\varphi _m^{(i)}H_m^{(kij)}J^{(j)},k=0,1,\mathrm{}N.`$ (14)
Below we prove the collinear factorization of the $`O(\alpha _s^{N+1})`$ diagrams $`G^{(N+1)}`$ into the convolution of the $`O(\alpha _s^N)`$ diagrams $`G^{(N)}`$ with the $`O(\alpha _s)`$ jet function $`J^{(1)}`$. Look for the gluon in a complete set of $`O(\alpha _s^{N+1})`$ diagrams $`G^{(N+1)}`$, one of whose ends attaches the lower most vertex on the internal quark line. Such a gluon exists, because $`G^{(N+1)}`$ are the finite-order diagrams. Let $`\alpha `$ denote the lower most vertex, and $`\beta `$ denote the attachments of the other end of the identified gluon inside the rest of the diagrams. There are two types of collinear configurations associated with this gluon, depending on whether the vertex $`\beta `$ is located on an internal line with a momentum along $`P_2`$. The fermion propagator adjacent to the vertex $`\alpha `$ is proportional to $`\mathit{}_2`$ in the collinear region with the loop momentum $`l`$ parallel to $`P_2`$. If $`\beta `$ is not located on a collinear line along $`P_2`$, the component $`\gamma ^{}`$ in $`\gamma ^\alpha `$ and the plus component of the vertex $`\beta `$ give the leading contribution. If $`\beta `$ is located on a collinear line along $`P_2`$, $`\beta `$ can not be plus, and both $`\alpha `$ and $`\beta `$ label the transverse components. This configuration is the same as of the self-energy correction to an on-shell particle.

According to the above classification, we decompose the tensor $`g_{\alpha \beta }`$ appearing in the propagator of the identified gluon into

$`g_{\alpha \beta }={\displaystyle \frac{n_{+\alpha }l_\beta }{n_+l}}\delta _\alpha \delta _\beta +\left(g_{\alpha \beta }{\displaystyle \frac{n_{+\alpha }l_\beta }{n_+l}}+\delta _\alpha \delta _\beta \right).`$ (15)
The first term on the right-hand side of Eq. (15) extracts the first type of collinear divergences, since the light-like vector $`n_{+\alpha }`$ selects the minus component of $`\gamma ^\alpha `$, and $`l_\beta `$ in the collinear region selects the plus component of the vertex $`\beta `$ . The second term extracts the second type of collinear divergences. The last term does not contribute a collinear divergence. We shall concentrate on the factorization corresponding to the first term, and the factorization corresponding to the second term can be achieved simply by requiring gauge invariance .

The identified collinear gluon with $`\alpha =`$ and $`\beta =+`$ does not attach the internal quark line directly, which carries a momentum along $`P_2`$ at the end point. That is, those diagrams with Fig. 1(e) as the $`O(\alpha _s)`$ subdiagram are excluded from the set of $`G^{(N+1)}`$ as discussing the first type of collinear configurations. Applying the Fierz transformation in Eq. (5) to break the fermion flow, we have the physical amplitude, in which the two on-shell quarks and the on-shell gluon carry the momenta $`\xi P_2`$, $`\xi P_2l`$ and $`l`$, respectively. The fraction $`\xi `$ reflects that the momentum flowing through the internal quark is almost parallel to $`P_2`$. Figure 2(a), describing the Ward identity, contains a complete set of contractions of $`l_\beta `$ represented by arrows, since the second diagram has been added back. The cuts on the internal quark lines denote the insertion of the Fierz identity.

The second diagram in Fig. 2(a) gives

$`l_\beta \xi \mathit{}_2\gamma ^\beta {\displaystyle \frac{1}{\xi \mathit{}_2\mathit{}}}=\xi \mathit{}_2(\mathit{}\xi \mathit{}_2+\xi \mathit{}_2){\displaystyle \frac{1}{\xi \mathit{}_2\mathit{}}}=\xi \mathit{}_2,`$ (16)
where the factor $`\xi \mathit{}_2`$ at the beginning of the above expression comes from the internal quark propagator adjacent to the photon vertex, and the term $`\xi \mathit{}_2`$ at the end leads to the $`O(\alpha _s^N)`$ diagrams. The diagrams $`G_{}^{(N+1)}`$ corresponding to the first term in Eq. (15) are then factorized according to Fig. 2(b). The factor $`n_{+\alpha }/n_+l`$ is just the Feynman rule associated with the Wilson line along the vector $`n_+`$ in Eq. (10). We arrive at the convolution for the first type of collinear configurations,

$`G_{}^{(N+1)}`$ $``$ $`G^{(N)}J_{}^{(1)},`$ (17)
with $`J_{}^{(1)}`$ given in Eq. (7).

The above procedures are applicable to the $`O(\alpha _s^{j+1})`$ jet function $`J^{(j+1)}`$. We identify the gluon in a complete set of $`O(\alpha _s^{j+1})`$ diagrams $`J^{(j+1)}`$, one of whose ends attaches the lower most vertex $`\alpha `$ on the internal quark line. The other end attaches the vertex $`\beta `$ inside the rest of the diagrams. For the first term on the right-hand side of Eq. (15), there exists a Ward identity similar to Fig. 2(a). It is then trivial to obtain the factorization of the jet function corresponding to the first type of collinear configurations,

$`J_{}^{(j+1)}`$ $``$ $`J^{(j)}J_{}^{(1)}.`$ (18)
The rest part of the proof is exactly the same as of the factorization theorem for the semileptonic decay $`B\pi l\overline{\nu }`$ in Sec. IV B of , with the pion distribution amplitudes replaced by the jet function, both of which are of the collinear origin. Employing Eqs. (14), (17) and (18) and following the steps of factorizing the decay $`B\pi l\overline{\nu }`$, we derive

$`G^{(N+1)}={\displaystyle \underset{m=+,}{}}{\displaystyle \underset{i=0}{\overset{N+1}{}}}{\displaystyle \underset{j=0}{\overset{N+1i}{}}}\varphi _m^{(i)}H_m^{(N+1ij)}J^{(j)},`$ (19)
with the $`O(\alpha _s^{N+1})`$ hard amplitude $`H_m^{(N+1)}`$ being infrared finite. Equation (19) indicates that all the soft and collinear divergences at the end point in the radiative decay $`B\gamma l\overline{\nu }`$ can be factorized into the $`B`$ meson distribution amplitudes and the jet function, respectively, order by order. We complete the proof of the factorization theorem, which is graphically described in Fig. 2(c), and construct the jet function defined in Eq. (10).

## IV JET FUNCTION

We now perform threshold resummation of the double logarithms $`\alpha _s\mathrm{ln}^2x`$ in covariant gauge $`A=0`$, which have been collected into the jet function to all orders in the previous section. Threshold resummation for inclusive QCD processes has been studied intensively . Here we shall adopt the framework developed in , which has been shown to lead to the same results as in . First, we replace the vector $`n_+`$ by $`n`$, which contains a plus component $`n^+>0`$ and a (small) minus component. This replacement, regularizing the collinear pole, extracts the double logarithm as stated in Sec. II. The jet functions in terms of the dimensionless vectors $`n_+`$ and $`n`$ can be regarded as being defined in the different factorization schemes. The definition then involves three variable vectors: the Wilson line direction $`n`$, the large momentum $`P_2`$, and the spectator momentum $`k`$. We argue that the jet function depends on the Lorentz invariants, $`nk`$, $`nP_2`$ and $`P_2k`$, since the other invariants $`P_2^2`$ and $`k^2`$ vanish, and $`n^2`$ will be fixed below. The third invariant $`P_2k`$ is not independent, because it can be rewritten as $`nP_2nk`$ with $`n^2`$ being a constant. We further argue that the scale invariance in $`n`$, as indicated by the Feynman rule associated with the eikonal line along $`n`$, implies that the jet function must depend on $`k`$ through the ratio $`nk/nP_2`$.

The next step is to derive the evolution of the jet function in $`x`$, i.e., in $`k^+=xP_1^+`$ by considering the derivative,

$`k^+{\displaystyle \frac{dJ}{dk^+}}={\displaystyle \frac{nk}{P_2k}}P_2^\alpha {\displaystyle \frac{dJ}{dn^\alpha }},`$ (20)
where we have applied the chain rule to relate the derivatives with respect to $`k`$ and to $`n`$. The differentiation $`d/dn^\alpha `$ operates on the eikonal line along $`n`$, giving

$`{\displaystyle \frac{nk}{P_2k}}P_2^\alpha {\displaystyle \frac{d}{dn^\alpha }}{\displaystyle \frac{n_\mu }{nl}}={\displaystyle \frac{\widehat{n}_\mu }{nl}},`$ (21)
with the special vertex,

$`\widehat{n}_\mu ={\displaystyle \frac{nk}{P_2k}}{\displaystyle \frac{P_2l}{nl}}n_\mu .`$ (22)
Another term in Eq. (22), proportional to $`P_2`$, is negligible, as the special vertex attaches inside of the jet function, which is dominated by momenta parallel to $`P_2`$.

The loop momentum $`l`$ flowing through the special vertex does not generate a collinear divergence due to vanishing of the numerator $`P_2l`$ in this region. It is easy to confirm that the ultraviolet region of $`l`$ does not produce $`\mathrm{ln}x`$ either. Therefore, we concentrate on the factorization of the soft gluon emitted from the special vertex, which can be achieved by applying the eikonal approximation to internal quark propagators, leading to $`n_\nu /n_{}l`$. Following the reasoning in , the derivative of the jet function is written as

$`x{\displaystyle \frac{dJ(x)}{dx}}=ig^2C_F{\displaystyle \frac{d^4l}{(2\pi )^4}\frac{\widehat{n}_\mu }{nl}\frac{g^{\mu \nu }}{l^2}\frac{n_\nu }{n_{}l}J(xl^+/P_1^+)},`$ (23)
where the argument of $`J`$ in the integral arises from the invariant mass of the internal quark, $`(P_2k+l)^22(xl^+/P_1^+)P_1P_2`$. Performing the integration over $`l^{}`$ and $`l_T`$, we derive the evolution equation,

$`x{\displaystyle \frac{dJ(x)}{dx}}={\displaystyle \frac{\alpha _s}{2\pi }}C_F{\displaystyle _x^1}{\displaystyle \frac{d\xi }{\xi x}}J(\xi ),`$ (24)
with the integration variable $`\xi `$ corresponding to the plus component $`l^+`$.

Employing the variable change $`\xi =x/z`$, Eq. (24) becomes

$`x{\displaystyle \frac{dJ(x)}{dx}}=J^{(0)}(x)+{\displaystyle \frac{\alpha _s}{2\pi }}C_F{\displaystyle _x^1}{\displaystyle \frac{dz}{z}}J(x/z),`$ (25)
where the nonperturbative inhomogenious term $`J^{(0)}(x)`$ absorbs the soft pole from $`\xi x`$ ($`z1`$). After removing this pole, the integral in the second term is infrared finite. It is interesting to observe that Eq. (25) is similar to the evolution equation for unintegrated parton distribution functions involved in inclusive QCD processes , which resums the same double logarithm $`\alpha _s\mathrm{ln}^2x`$.

To solve Eq. (24), we perform the Mellin transformation from the momentum fraction ($`x`$) space to the moment ($`N`$) space:

$`\stackrel{~}{J}(N){\displaystyle _0^1}𝑑x(1x)^{N1}J(x).`$ (26)
The small $`x`$ region then corresponds to the large $`N`$ region. The left-hand side of Eq. (24) leads to

$`{\displaystyle _0^1}𝑑x(1x)^{N1}x{\displaystyle \frac{dJ(x)}{dx}}=N\stackrel{~}{J}(N)+(N1)\stackrel{~}{J}(N1){\displaystyle \frac{d}{dN}}\left[N\stackrel{~}{J}(N)\right],`$ (27)
which is reasonable under the approximation $`dN\mathrm{\Delta }N=1`$. The right-hand side is given by

$`{\displaystyle \frac{\alpha _s}{2\pi }}C_F{\displaystyle _0^1}𝑑\xi (1\xi )^{N1}J(\xi ){\displaystyle _0^1}{\displaystyle \frac{dw}{1w}}\left({\displaystyle \frac{1\xi w}{1\xi }}\right)^{N1}{\displaystyle \frac{\alpha _s}{2\pi }}C_F\mathrm{ln}N\stackrel{~}{J}(N),`$ (28)
where we have exchanged the sequence of the integrations over $`x`$ and over $`\xi `$, made the variable change $`x=\xi w`$, and kept only the $`N`$-dependent term. The trick for extracting the $`N`$ dependence from the integral over $`w`$ in Eq. (28) is to consider its derivative with respect to $`N`$, and to insert the large $`N`$ approximation for the $`\delta `$-function, $`lim_N\mathrm{}N(1y)^{N1}=\delta (y)`$ with $`y=\xi (1w)/(1\xi w)`$.

Equating Eqs. (27) and (28), we have the evolution equation in the moment space,

$`N{\displaystyle \frac{d\stackrel{~}{J}(N)}{dN}}=\stackrel{~}{J}(N){\displaystyle \frac{\alpha _s}{2\pi }}C_F\mathrm{ln}N\stackrel{~}{J}(N),`$ (29)
whose solution is the Sudakov factor,

$`\stackrel{~}{J}(N)={\displaystyle \frac{1}{N}}\mathrm{exp}\left({\displaystyle \frac{1}{4}}\gamma _K\mathrm{ln}^2N\right),`$ (30)
with the anomalous dimension $`\gamma _K=\alpha _sC_F/\pi `$. The factor $`1/N`$ can be regarded as the Mellin transformation of the initial condition $`J^{(0)}=1`$. Comparing Eq. (30) with the resummation for an inclusive jet , the resummation for an exclusive jet is suppressed by a factor $`1/N`$, and the corresponding anomalous dimension is down by a factor 2.

By means of the inverse Mellin transformation, the jet function in the momentum fraction space is written as

$`J(x)={\displaystyle _{ci\mathrm{}}^{c+i\mathrm{}}}{\displaystyle \frac{dN}{2\pi i}}(1x)^N\stackrel{~}{J}(N),`$ (31)
with $`c`$ being an arbitrary real constant larger than all the real parts of the poles of the integrand. Since $`\stackrel{~}{J}(N)`$ contains a branch cut along the negative real axis on the complex $`N`$ plane, Eq. (31) reduces to

$`J(x)=\mathrm{exp}\left({\displaystyle \frac{\pi }{4}}\alpha _sC_F\right){\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dt}{\pi }}(1x)^{\mathrm{exp}(t)}\mathrm{sin}\left({\displaystyle \frac{1}{2}}\alpha _sC_Ft\right)\mathrm{exp}\left({\displaystyle \frac{\alpha _s}{4\pi }}C_Ft^2\right),`$ (32)
where the variable change $`N=\mathrm{exp}(t+i\pi )`$ ($`N=\mathrm{exp}(ti\pi )`$) has been adopted for the piece of contour above (below) the branch cut as shown in Fig. 3. Note that the above expression holds only for $`\alpha _s>0`$. For $`\alpha _s=0`$, the residue associated with the $`N=0`$ pole should be included. It is trivial to check that $`J(x)`$ is normalized to unity, $`J(x)𝑑x=\stackrel{~}{J}(1)=1`$ .

Obviously, Eq. (32) vanishes at $`x0`$, because the integrand is an odd function in $`t`$, and at $`x1`$ due to the factor $`(1x)^{\mathrm{exp}(t)}`$. The latter property is the consequence of the extrapolation of the Sudakov factor to the low $`N`$ region. Moreover, Eq. (32) provides suppression near the end point $`x0`$, which is stronger than any power of $`x`$. This is understood from vanishing of all the derivatives of Eq. (32) with respect to $`x`$ at $`x0`$. For example, the first derivative gives

$`{\displaystyle \frac{d}{dx}}J(x)|_{x0}`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{\pi }{4}}\alpha _sC_F+{\displaystyle \frac{\pi }{\alpha _sC_F}}\right){\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dt}{\pi }}e^t\mathrm{sin}\left({\displaystyle \frac{1}{2}}\alpha _sC_Ft\right)\mathrm{exp}\left({\displaystyle \frac{\alpha _s}{4\pi }}C_Ft^2\right),`$ (33)
where the variable change $`tt+2\pi /(\alpha _sC_F)`$ has been made. The integrand in Eq. (33) is also an odd function in $`t`$, and the integral diminishes. To the accuracy of the next-to-leading logarithms, the running of the coupling constant $`\alpha _s`$ should be taken into account, and Eq. (32) will be modified. However, the above features remain.

## V SEMILEPTONIC DECAY

In this section we extend the above formalism to the semileptonic decay $`B\pi l\overline{\nu }`$ in the fast recoil region of the pion. The $`B`$ meson momentum $`P_1`$ is the same as in the decay $`B\gamma l\overline{\nu }`$, and the pion momentum $`P_2`$ is the same as the photon momentum. Leading-twist factorization theorem for the $`B\pi `$ form factor $`f_+(q^2)`$ in the standard definition has been proved in ,

$`f_+(q^2)={\displaystyle \underset{m=+,}{}}\varphi _m(x_1)H_m(x_1,x_2,\eta )\varphi _\pi (x_2),`$ (34)
which holds in the region with $`x_1O(\overline{\mathrm{\Lambda }}/M_B)`$ and with $`x_2O(1)`$.

Since Fig. 4(a), proportional to $`1/(x_1x_2^2)`$, is more singular at small $`x_2`$, we consider the end-point region with $`x_2O(\overline{\mathrm{\Lambda }}/M_B)`$, where the internal $`b`$ quark propagator scales like $`1/(\overline{\mathrm{\Lambda }}M_B)`$. Part of $`O(\alpha _s)`$ corrections to Fig. 4(a) are shown in Figs. 4(c)-(e), among which Fig. 4(c) generates the double logarithm $`\alpha _s\mathrm{ln}^2x_2`$ from the collinear region with the loop momentum parallel to $`P_2`$. Figures 4(d) and 4(e) involve only the single logarithm. We emphasize that the double logarithm discussed here for Fig. 4(a) differs from that in , which concerns $`\alpha _s\mathrm{ln}^2x_1`$. As stated above, Fig. 4(a) is less singular at $`x_10`$, and $`\alpha _s\mathrm{ln}^2x_1`$ is not relevant. The double logarithm in Fig. 4(c) is simply extracted by eikonalizing the light quark line, which flows into the pion, and the $`b`$ quark line according to Eq. (3):

$`J_{}^{(1)}`$ $`=`$ $`ig^2C_F{\displaystyle \frac{d^4l}{(2\pi )^4}\frac{n_+n_{}}{n_+(lx_2P_2)n_{}ll^2}},`$ (35)
$`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}C_F\mathrm{ln}^2x_2+\mathrm{}.`$ (36)
Similarly, the collinear pole has been regularized by including a small plus component for $`n_{}`$. It is observed that the double logarithm in Eq. (36) is the same as in Eq. (7).

The all-order factorization of the jet function from the decay $`B\gamma l\overline{\nu }`$ can be proved following the procedures in Sec. III. It is also easy to show that this jet function obeys the evolution equation in Eq. (24) by identifying the correspondence of $`x_2P_2`$ and $`n_+`$ to $`k`$ and $`P_2`$ in Sec. IV, respectively. Hence, the threshold resummation leads to a result the same as Eq. (32). That is, the Sudakov factor is universal.

The analysis for Fig. 4(b) is similar to that for the decay $`B\gamma l\overline{\nu }`$. In the end-point region with $`x_1O(\overline{\mathrm{\Lambda }}^2/M_B^2)`$, additional collinear divergences associated with the internal light quark are produced in Figs. 4(f)-4(h). Figure 4(f) gives the double logarithm $`\alpha _s\mathrm{ln}^2x_1`$, whose factorization is the same as of Fig. 1(c). Note that Fig. 5 has been considered as the dominant source of the double logarithm $`\alpha _s\mathrm{ln}^2x_1`$ associated with Fig. 4(b) . It will be demonstrated explicitly in the Appendix A that Fig. 5 in fact contains only the single logarithm $`\alpha _s\mathrm{ln}x_1`$. The result in  might be attributed to an inappropriate collinear approximation. Therefore, the universality of the double-logarithm resummation indeed holds.

At last, we investigate the resummation effect on the $`B\pi `$ form factor $`f_+(q^2)`$, whose factorization formula is given by

$`f_+(q^2)`$ $`=`$ $`{\displaystyle \frac{\pi \alpha _sC_Ff_Bf_\pi }{\eta ^2M_B^2N_c}}{\displaystyle 𝑑x_1𝑑x_2\varphi _B(x_1)\left[\frac{J(x_2)}{x_1x_2^2}(1+x_2\eta )\frac{J(x_1)}{x_1x_2}(1\eta )\right]\varphi _\pi (x_2)},`$ (37)
with $`N_c=3`$ being the number of colors. We employ the models ,

$`\varphi _\pi (x)`$ $`=`$ $`6x(1x)\{1+0.66[5(12x)^21]\},`$ (38)
$`\varphi _B(x)`$ $``$ $`{\displaystyle \frac{1}{2}}[\varphi _+(x)\varphi _{}(x)]=N_B\sqrt{x(1x)}\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{xM_B}{\omega _B}}\right)^2\right],`$ (39)
with the shape parameter $`\omega _B`$ and the normalization constant $`N_B`$ determined by $`\varphi _B(x)𝑑x=1`$. Another $`B`$ meson distribution amplitude $`\overline{\varphi }_B(x)=[\varphi _+(x)\varphi _{}(x)]/2`$ is power-suppressed, since its first nonvanishing moment starts with $`O(1/M_B)`$ . It has been confirmed, contrary to the conclusion in , that the contribution from $`\overline{\varphi }_B`$ is less important than that from $`\varphi _B`$ . For the purpose of this work, it is sufficient to work on a single $`B`$ meson distribution amplitude. If $`J(x)`$ is excluded, the first term in Eq. (37) is logarithmically divergent. With the threshold resummation, the form factor is calculable without introducing any infrared cutoffs .

Choosing $`\alpha _s=0.4`$, $`M_B=5.28`$ GeV, and the decay constants $`f_B=190`$ MeV and $`f_\pi =130`$ MeV, we obtain $`f_+(q^2)`$ for $`\omega _B=0.3(0.4)`$ GeV,

$`f_+(0)=0.15(0.12),f_+(2\mathrm{GeV}^2)=0.18(0.13),f_+(4\mathrm{GeV}^2)=0.21(0.15).`$ (40)
The variation of $`f_+`$ for $`\alpha _s=0.3`$-0.5 is less than 10%, because the change of $`J`$ is compensated by that of $`\alpha _s`$ in the overall coefficient in Eq. (37). The difference between the above values and the expected one $`f_+(0)0.3`$ can be resolved by taking into account higher-twsit contributions . Including two-parton twist-3 pion distribution amplitudes, which are finite at the end point , the singularities in the $`B\pi `$ form factor become linear. In this case threshold resummation is even more crucial.

## VI CONCLUSION

In this paper we have shown that the double logarithmic corrections $`\alpha _s\mathrm{ln}^2x`$ appear in exclusive $`B`$ meson decays. When the end-point region with a momentum fraction $`x0`$ is important, these double logarithms need to be organized to all orders in order to justify perturbative expansion. The double logarithms, associated with internal particles which almost go on mass shell, are of the collinear nature. The PQCD factorization theorem for exclusive $`B`$ meson decays at the end point then demands the introduction of a quark jet function into decay amplitudes. The factorization of the jet function from the $`B\gamma l\overline{\nu }`$ mode has been proved rigorously. The proof for the factorization from other modes, such as the semileptonic decay $`B\pi l\overline{\nu }`$, is similar. It has been shown that the jet function is defined as a matrix element of a quark field attached by a Wilson line, based on which threshold resummaiton can be performed. It is interesting to achieve the above proof in the framework of the soft-collinear effective theory .

We have derived the evolution equation for the jet function, which is solvable in the large $`N`$ limit in the Mellin space. The solution is a universal, i.e., process-independent, Sudakov factor. The qualitative behavior of this Sudakov factor has analyzed and found to decrease quickly at $`x0`$. It has been demonstrated that Sudakov suppression is strong enough to smear the end-point singularity in the $`B\pi `$ form factor and leads to reasonable predictions. We conclude that in a self-consistent PQCD analysis of the heavy-to-light transition form factors, the end-point singularities do not exist. In a future work we shall extend the formalism developed here to more complicated nonleptonic $`B`$ meson decays. Threshold resummation for various topologies of decay amplitudes, such as annihilation and nonfactorizable ones, will be studied.

I thank M. Beneke, Y.Y. Keum, T. Kurimoto, D. Pirjol, A.I. Sanda, and G. Sterman for useful discussions. The work was supported in part by the National Science Council of R.O.C. under Grant No. NSC-90-2112-M-001-077, by the National Center for Theoretical Sciences of R.O.C., by Grant-in Aid for Special Project Research (Physics of CP violation), and by Grant-in Aid for Scientific Exchange from the Ministry of Education, Science and Culture of Japan.

## A TRIPLE-GLUON CORRECTION

In this Appendix we show that Fig. 5 produces only the single logarithm $`\alpha _s\mathrm{ln}x_1`$. The loop integral is written as

$`I`$ $`=`$ $`g^4{\displaystyle }{\displaystyle \frac{d^4l}{(2\pi )^4}}tr[\gamma _\delta {\displaystyle \frac{\gamma _5\mathit{}_2}{4N_c}}\gamma _\beta {\displaystyle \frac{\mathit{}_2\mathit{}_1+\mathit{}}{(P_2k_1+l)^2}}\gamma ^\mu `$ (A3)
$`\times {\displaystyle \frac{\mathit{}_1\mathit{}_1+\mathit{}+M_B}{(P_1k_1+l)^2M_B^2}}\gamma _\alpha {\displaystyle \frac{(\mathit{}_1+M_B)\gamma _5}{4N_c}}]`$
$`\times {\displaystyle \frac{\mathrm{\Gamma }^{\alpha \beta \delta }f^{abc}tr(T^aT^bT^c)}{(k_1x_2P_2l)^2l^2(x_2P_2k_1)^2}},`$
with the triple-gluon vertex,

$`\mathrm{\Gamma }^{\alpha \beta \delta }=(2x_2P_22k_1+l)^\alpha g^{\beta \delta }+(k_1x_2P_22l)^\delta g^{\alpha \beta }+(lx_2P_2+k_1)^\beta g^{\alpha \delta },`$ (A4)
the anti-symmetric tensor $`f^{abc}`$, and the color matrices $`T^{a,b,c}`$.

We concentrate on the frist term in Eq. (A4), since the calculation of the second and third terms is similar. If the loop momentum $`l`$ in the quark propogator $`1/(P_2k_1+l)^2`$ is neglected, Eq. (A3) gives the double logarithm $`\alpha _s\mathrm{ln}^2x_1`$ as obtained in . However, this approximation is not appropriate for the collection of the collinear divergence, for which $`l`$ is of the same order as $`P_2`$. Keeping the $`l`$ dependence, and using the identity $`f^{abc}tr(T^aT^bT^c)=6i`$, we have

$`I_1`$ $`=`$ $`{\displaystyle \frac{ig^4}{4N_c}}{\displaystyle \frac{tr(\mathit{}_2\mathit{}_1\gamma ^\mu \mathit{}_1)}{(x_2P_2k_1)^2}}`$ (A6)
$`\times {\displaystyle }{\displaystyle \frac{d^4l}{(2\pi )^4}}{\displaystyle \frac{v(2x_2P_22k_1+l)}{vl(P_2k_1+l)^2(k_1x_2P_2l)^2l^2}}.`$
The contribution proportional to $`M_B^2`$ in Eq. (A3) is suppressed by a power of $`x_1`$ compared to Eq. (A6). For a similar reason, the term $`vk_1`$ in Eq. (A6) is also negligible. The term $`vl`$ in the numerator does not lead to a collinear divergence, because the poles of the corresponding denominator $`(P_2k_1+l)^2(k_1x_2P_2l)^2l^2`$ in the $`l^{}`$ plane are not of the pinched type.

Doing the contour integral with the residue $`l^{}=l^+iϵ`$ (for $`l^+<0`$), Eq. (A6) becomes

$`I_1`$ $`=`$ $`H_{4b}^{(0)\mu }{\displaystyle \frac{4g^2}{(2\pi )^3C_F}}{\displaystyle _{\mathrm{}}^0}𝑑l^+{\displaystyle _{\mathrm{}}^{\mathrm{}}}d^2l_T`$ (A8)
$`\times {\displaystyle \frac{x_2P_2^{}(P_2k_1)^2}{(l_T^2+2l^{+2}2P_2^{}l^++2P_2^{}k_1^+)(l_T^2+2l^{+2}2x_2P_2^{}l^++2x_2P_2^{}k_1^+)(l_T^2+2l^{+2})}},`$
where

$`H_{4b}^{(0)\mu }={\displaystyle \frac{g^2C_Ftr(\mathit{}_2\mathit{}_1\gamma ^\mu \mathit{}_1)}{8N_c(x_2P_2k_1)^2(P_2k_1)^2}},`$ (A9)
is the lowest-order amplitude from Fig. 4(b). The integration over $`l_T`$ leads to

$`I_1`$ $`=`$ $`H_{4b}^{(0)\mu }{\displaystyle \frac{\alpha _s}{4\pi }}C_Ax_2k_1^+{\displaystyle _0^{\mathrm{}}}dl^+[{\displaystyle \frac{\mathrm{ln}(l^{+2}/M_B^2+\eta l^+/M_B+x_1\eta )}{(1x_2)(l^++k_1^+)^2}}`$ (A11)
$`{\displaystyle \frac{\mathrm{ln}(l^{+2}/M_B^2+x_2\eta l^+/M_B+x_1x_2\eta )}{x_2(1x_2)(l^++k_1^+)^2}}+{\displaystyle \frac{\mathrm{ln}(l^{+2}/M_B^2)}{x_2(l^++k_1^+)^2}}],`$
with the color factor $`C_A=3`$.

The terms $`l^{+2}/M_B^2`$ in the first two logarithms, giving contributions suppressed by a power of $`x_1`$ compared to those from $`l^+/M_B`$, can be dropped. Each term in Eq. (A11) then produces only the single logarithms $`\alpha _s\mathrm{ln}x_1`$:

$`I_1`$ $`=`$ $`H_{4b}^{(0)\mu }{\displaystyle \frac{\alpha _s}{4\pi }}C_Ax_2\left[{\displaystyle \frac{\mathrm{ln}x_1}{1x_2}}{\displaystyle \frac{\mathrm{ln}x_1}{x_2(1x_2)}}+2{\displaystyle \frac{\mathrm{ln}x_1}{x_2}}\right],`$ (A12)
$`=`$ $`H_{4b}^{(0)\mu }{\displaystyle \frac{\alpha _s}{4\pi }}C_A\mathrm{ln}x_1.`$ (A13)
In the above expression the constant terms independent of $`\mathrm{ln}x_1`$ are not displayed. It can be shown in a similar way that the other two terms in Eq. (A4) involve only the single logarithms.

Figure Captions

Fig. 1. (a) and (b) Lowest-order diagrams for the decay $`B\gamma l\overline{\nu }`$. (c)-(h) $`O(\alpha _s)`$ corrections to Fig. 1(a).

Fig. 2. (a) The Ward identity. (b) Factorization of the $`O(\alpha _s)`$ diagrams as a consequence of (a). (c) Factorization of the $`B\gamma l\overline{\nu }`$ decay amplitude in the end-point region.

Fig. 3. Contour for the inverse Mellin transformation of the jet function.

Fig. 4. (a) and (b) Lowest-order diagrams for the decay $`B\pi l\overline{\nu }`$. (c)-(e) \[(f)-(h)\] Part of $`O(\alpha _s)`$ corrections to Fig. 4(a) \[Fig. 4(b)\].

Fig. 5. Triple-gluon correction to Fig. 4(b).