# Evolution of Compact-Binary Populations

in Globular Clusters: A Boltzmann Study

I. The Continuous Limit

###### Abstract

We explore a Boltzmann scheme for studying the evolution of compact binary populations of globular clusters. We include processes of compact-binary formation by tidal capture and exchange encounters, binary destruction by dissociation and other mechanisms, and binary hardening by encounters, gravitational radiation and magnetic braking, as also the orbital evolution during mass transfer, following Roche lobe contact. For the encounter processes which are stochastic in nature, we study the probabilistic, continuous limit in this introductory work, deferring the specific handling of the stochastic terms to the next step. We focus on the evolution of (a) the number of X-ray sources in globular clusters, and (b) the orbital-period distribution of the X-ray binaries, as a result of the above processes. We investigate the dependence of on two essential cluster properties, namely, the star-star and star-binary encounter-rate parameters and , which we call Verbunt parameters. We compare our model results with observation, showing that the model values of and their expected scaling with the Verbunt parameters are in good agreement with results from recent X-ray observations of Galactic globular clusters, encouraging us to build more detailed models.

^{†}

^{†}slugcomment: To appear in Ap. J.

## 1 Introduction

In this era of high-resolution X-ray observations with *Chandra*
and *XMM-Newton*, studies of compact binaries in globular clusters
have reached an unprecedented level of richness and detail. The numbers
of compact X-ray binaries detected in Galactic globular clusters with
high central densities are now becoming large enough that diagnostic
correlations with essential cluster parameters, such as the two-body
encounter rate , can be performed (Pooley et al., 2003) at a high
level of statistical significance. The results of such observational
studies are naturally to be compared with those obtained from
theoretical modeling of binary dynamics in globular clusters, which
has had a long history, from the pioneering semi-analytic work of the
1970s (Heggie, 1975), to the more detailed numerical scattering experiments
of the 1980s (Hut & Bahcall, 1983), leading to the wealth of detailed numerical
work of the early- to mid-1990s (Makino & Aarseth, 1992; Heggie & Hut, 2003) using a variety of
techniques including Fokker-Planck and Monte Carlo approaches, as also
N-body simulations, and finally to the extensive N-body simulations in
the latter half of the 1990s using special-purpose computers with
ultrahigh speeds (Makino & Taiji, 1998; Hut, 2001).

The range of problems studied by the above modeling has also been
extensive. From the study and classification of individual scattering
events to the construction of comprehensive fitting formulas for the
cross-sections of such events (Hut & Bahcall, 1983; Heggie, Hut & McMillan, 1996), from the development
of Fokker-Planck codes to the use of Monte Carlo methods for following
binary distributions in globular clusters (Gao et.al., 1991; Hut et.al., 1992), and
from tracking the fate of a relatively modest population of test
binaries against a fixed stellar background to being able to tackle
similar projects for much larger binary populations with the aid of the
above special-purpose machines (Hut et.al., 1992; Makino, 1996), efforts along
various lines of approach have shed light on the overall phenomenon of
binary dynamics and evolution in globular clusters from various angles.
For example, evolutions of the distributions of both *external* and
*internal* binding energies of the binaries under stellar encounters
have been studied by several authors, the emphasis usually being on the
former, and final results on the external binding energy being expressed
almost universally in terms of their radial positions inside the
cluster, which provides an equivalent description
(Hut et.al., 1992; Sigurdsson & Phinney, 1993, 1995).

In this series of papers, we introduce an alternative method of studying the
evolution of compact-binary populations in globular clusters, wherein we
use a Boltzmann description to follow the time-evolution of such populations,
subject to both (a) those processes which determine compact-binary evolution
in isolation (i.e., outside globular clusters, or, in the “field” of the
host galaxy, so to speak), e.g., angular momentum loss by gravitational
radiation and magnetic braking, as also orbital evolution due to mass
transfer, and, (b) those processes which arise from encounters of compact
binaries with the dense stellar background in globular clusters, e.g.,
collisional hardening (Heggie, 1975; Shull, 1979; Banerjee & Ghosh, 2006), binary formation through
tidal capture and exchange processes, and binary destruction. We treat all
of the above processes simultaneously through a Boltzmann formalism, the aim
being to see their combined effect on the compact-binary population as a
whole, in particular on the evolution of (a) the total number of X-ray
binaries as the formation and destruction processes continue to operate,
and, (b) the orbital-period distribution of the population. We stress at
the outset that ours is not a Fokker-Planck description but the original
Boltzmann one, which in principle is capable of handling both the
*combined* small effects of a large number of frequent, weak,
distant encounters *and* the *individual* large effects of a
small number of rare, strong, close encounters. In our approach, both of
the above two types of effects are taken into account through
cross-sections for the relevant processes, as determined from extensive
previous work on numerical experiments with two-body and three-body
encounters (Heggie, Hut & McMillan, 1996; Portegies Zwart et.al., 1997b). As these processes are inherently
stochastic, a natural question that arises is how they are to be handled
simultaneously with those which govern the fate of isolated compact
binaries, and which are inherently continuous. It is essential to
appreciate the importance of this question, since a simultaneous action
of the above continuous and stochastic processes is precisely what
operates on binaries in globular clusters, and so produces the observed
properties of compact-binary populations in them.

Our answer to the above question is a step-by-step one, as follows. As the first step, in this first paper of the series (henceforth Paper I), we explore the continuous limit of the above stochastic processes, wherein the probability or cross-section of a particular such process happening with a given set of input and output variables is treated as a continuous function of these variables. This is, of course, a simplification, but it serves as a clarification of the average, long-term trends expected in the evolution of the binary population. In the next step, in the second paper of the series (henceforth Paper II), we treat the stochastic processes as stochastic terms in the Boltzmann equation with cross-sections as given in Paper I, with the aid of relatively recently-developed methods for solving stochastic partial differential equations. The resulting evolutionary trends show stochastic behavior, as expected, with fluctuations that vary from one particular “realization” of the essential processes to another. However, the average trends follow the continuous limit computed in Paper I, which is as expected, and which shows the relevance of extracting this limit.

In Papers I and II, we model the stellar background provided by the globular cluster as a fixed background with given properties, as has been widely done in previous works (Hut et.al., 1992; Portegies Zwart et.al., 1997b; Sigurdsson & Phinney, 1993, 1995): this amounts to neglecting the back reaction of binary evolution on the background, which is reasonable if the main aim is an investigation of essential features of binary evolution, as was the case in the above previous works, as also in this work. However, the globular-cluster background does evolve slowly, passes through the core-collapse phase and possible gravothermal oscillations (Sugimoto & Bettwieser, 1983; Gao et.al., 1991), so that it would be interesting to be able to follow the effects of these on the evolution of the compact-binary population. We do this in the third paper of the series (henceforth Paper III), wherein we adopt previous results on time-evolution of globular-cluster properties, and study their effects on the evolution of compact-binary populations, again under the approximation of neglecting the back reaction of binary evolution on the globular-cluster background, as above and as appropriate for a first look.

In our study, we focus primarily on two aspects of the compact-binary populations of globular clusters. First, we study how the total number of X-ray binaries (henceforth XBs, which are mass-transferring compact binaries where the donor is a low-mass “normal” star, and the accretor is a degenerate star — a neutron star or a heavy white dwarf) in a cluster evolves as the stellar encounter processes proceed. Second, we also follow the evolution of the orbital-period () distribution of the pre-X-ray binaries (henceforth PXBs; also see below) and XBs, (or, equivalently, the distribution of their orbital radii ) within the framework of our model. However, we have adopted here only a very simple model of orbital evolution of individual binaries in order to assess the feasibility of our basic approach to globular-cluster environments, as detailed later. Consequently, while the -distribution found by us may be roughly applicable to cataclysmic variables (CVs) with white-dwarf accretors, it cannot be compared at this stage to that of low-mass X-ray binaries (LMXBs) with neutron-star accretors, without including the essential stellar evolutionary processes that occur during the PXB and XB phase. Thus, we record our computed -distribution here only as a preliminary indication of the results that emerge naturally from this line of study at this stage, to be improved upon later.

The basic motivation for our study comes from recent advances in X-ray
observations of globular clusters, as mentioned above: with sufficient
numbers of X-ray binaries detected in globular clusters, an understanding
of how is influenced by essential globular-cluster parameters is
becoming a central question. With the above goal in mind, we therefore
explicitly follow the evolution of binaries only in internal binding
energy (or binary period, or binary separation, which are equivalent
descriptions if the stellar masses are known) and time, but not of their
external binding energy (or position inside the globular cluster; see
above). We emphasize that we do not *neglect* changes in the latter
in any way, as they are automatically taken care of in the detailed dynamics
of encounters which are represented by the relevant cross-sections mentioned
above and elaborated on in the following sections. It is only that we do not
keep an explicit account of them, as we do not need them for our purposes.
In other words, we consider a bivariate binary distribution function
, which may be looked upon as the integral of the distribution
over all admissible values of , or
equivalently over all positions inside the globular cluster
(Hut et.al., 1992; Sigurdsson & Phinney, 1993, 1995). We also emphasize that, by doing so, we do
*not* implicitly assume any particular correlation, nor a lack thereof,
between and (Hut et.al., 1992): whatever correlations
result from the dynamics of the encounters will be automatically displayed
if we follow the evolution in or , which is not of interest to
us in this particular study.

Our first results from the above evolutionary scheme show that the total number of XBs expected in a globular cluster scales in a characteristic way with well-known globular cluster parameters and (which we call Verbunt parameters: see Sec. 2.1) whose qualitative nature is rather similar to that found in our earlier “toy” model (Banerjee & Ghosh, 2006), although some details are different. Basically, scales with — a measure of the dynamical formation rate of compact binaries, and, at a given , decreases with increasing at large values of — a measure of the rate of destruction of these binaries by dynamical processes. These expected theoretical trends with the Verbunt parameters compare very well with the observed trends in recent data, encouraging us to construct more detailed evolutionary schemes.

In Sec. 2, we detail our model of compact binary evolution in globular clusters, describing, in turn, our handling of globular clusters, binary formation, destruction, and hardening processes, our Boltzmann scheme for handling population-evolution, and our numerical method. In Sec. 3, we give our model results on (a) the expected number of X-ray binaries in globular clusters as a function of their Verbunt parameters, and (b) the evolution of compact-binary period distribution. In Sec. 4, we compare these model results with the current observational situation. Finally, we collect our conclusions and discuss future possibilities in Sec. 5.

## 2 Model of Compact Binary Evolution in Globular Clusters

We consider a binary population described by a number distribution , where is the binary separation, interacting with a fixed background of stars representing the core of a globular cluster of stellar density and core radius . We now describe various ingredients of our model and the evolutionary scheme.

### 2.1 Globular clusters

Globular cluster cores are described by an average stellar density ,
a velocity dispersion , and a core radius . In this work, we
consider star-star and star-binary encounters of various kinds, but neglect
binary-binary encounters. For characterizing the former two processes, two
encounter rates are defined and used widely (Verbunt, 2002, 2006).
The first is the two-body stellar encounter rate ,
which scales with , and occurs naturally in the rates of
two-body processes like tidal capture, stellar collisions and merger.
In fact, we can *define* it as

(1) |

for our purposes here. Note that the last scaling in the above equation holds only for virialized cores, where the scaling can be applied. In this work, we shall use this assumption where necessary, but with the caveat that some observed globular clusters have clearly not virialized yet.

The second is a measure of the rate of encounter between binaries and single
stars in the cluster, the rate normally used being the encounter rate
of a *single* binary with the stellar background, with the
understanding that the total rate of binary-single star encounter in the
cluster will be . We can *define* for our
purposes as we did in Banerjee & Ghosh (2006), namely,

(2) |

where the last scaling holds, again, only for virialized cores.

The importance of the above cluster parameters and in this
context has been extensively discussed by Verbunt (Verbunt, 2002, 2006), and
we shall call them *Verbunt parameters* here. Note that, for virialized
cores, we can invert Eqs. (1) and (2) to obtain
the scaling of the core density and radius with the Verbunt parameters as:

(3) |

It is most instructive to display the observed globular clusters in the
plane, which we do^{1}^{1}1Alternatively, the display can
be in the plane, as in Verbunt’s original work. We find the
cluster dynamics more transparent when shown directly in terms of the
Verbunt parameters. in Fig. 1. The point that immediately strikes
one in the figure is that the observed globular clusters seem to occur in a
preferred, diagonal, “allowed” band in the plane, along
which there is a strong, positive correlation between the two parameters.
We shall return to the significance of this elsewhere.

In Fig. 1, we also overplot the positions of those clusters in which significant numbers of X-ray sources have been detected, color-coding them according to the number of X-ray sources in each of them, as indicated. It is clear that these clusters are all in the upper parts of the above “allowed” band, which is entirely consistent with the widely-accepted modern idea that the dominant mechanisms for forming these compact XBs in globular clusters are dynamical, e.g., tidal capture, exchange encounters, and so on, since such mechanisms occur more efficiently at higher values of the Verbunt parameters and , corresponding to higher stellar densities in the cluster core. Note that the probability of destruction of binaries by dynamical processes also increases with increasing , as we shall see below, so that, at first sight, we might have expected the highest incidence of XBs in those clusters which have high and low . However, since and are strongly correlated positively, as above, we cannot have arbitrarily high and low for the same cluster. In reality, the highest number of XBs seem to occur, as Fig. 1 shows, in those clusters which have the highest values of and high, but not the highest, values of . We return to this point later in the paper, where we present our theoretical expectations for the scaling of the number of binary X-ray sources with the Verbunt parameters and on the basis of the evolutionary scheme explored here.

In modeling the globular cluster core as a static background in this work, we assume that, initially, a fraction of the stars is in primordial binaries, and that a fraction of the stellar population is compact, degenerate stars with the canonical mass (representing neutron stars and heavy white dwarfs). The rest of the stellar background (including the primordial binaries) is taken to consist of low-mass stars of the canonical mass , which is a reasonable estimate of the mean stellar mass of a mass-seggregated core (Portegies Zwart et.al., 1997a). Naturally, the compact binaries formed from these ingredients consist of a degenerate star of mass , and a low-mass companion of mass . While this is clearly an oversimplification which must be improved upon in subsequent work, it appears to be adequate for a first look, which is our purpose here.

### 2.2 A Boltzmann evolutionary scheme

We explore in this work a Boltzmann evolutionary scheme, wherein the evolution of the number of binaries per unit interval in the binary separation (we choose to work here with ; equivalent descriptions in terms of the binary period or the internal binding energy [see Sec. 1] are possible, of course) is described by

(4) |

Here, is the total derivative of bivariate : as explained in
Sec. 1, this is the result of an integration of a
general, multivariate binary distribution over the variables we do not
follow explicitly in this study, e.g., the external binding energy or,
equivalently, the position of the binary inside the globular cluster.
Further, is the total rate of binary formation per unit interval
in due to the various processes detailed below, and is the
total rate of binary destruction *per binary* per unit interval in
due to various processes, also detailed below. As our model stellar
background representing the cluster core is taken as static for Papers I
and II, the Verbunt parameters and are time-independent,
so that the formation and destruction rates and only depend on
and the stellar masses.

The above evolution equation can be re-written in the usual Boltzmann form

(5) |

where represents the total rate of shrinkage or
*hardening* of binaries (i.e., ) due to several effects, which
we introduced in Sec. 1, and which we elaborate on below. In the
absence of all processes of formation and destruction, ,
Eq. (5) becomes the usual collisionless Boltzmann equation

(6) |

representing a movement or “current” of binaries from larger to smaller values of due to hardening. Equation (6) as akin to a wave equation with a formal “phase velocity” of propagation. This analogy often proves useful for solving many problems, even with the more complicated formation and destruction terms present in Eq. (5). Note that, when is constant (or roughly so, which can happen under certain circumstances, as we shall see later), the elementary wave-equation analogy is quite exact, and solutions of the form should apply. We shall explore this point elsewhere.

Note further that the Boltzmann scheme outlined above does not have an explicit inclusion of the escape of those binaries from the globular cluster which receive a sufficiently large “kick”. In principle, we can include this by suitably generalizing the above destruction term . However, in this introductory study, this did not appear crucial, as the main population affected by this process is that of primordial binaries, whereas our main concern here is with dynamically-formed compact binaries. The latter are, generally speaking, already so hard at formation that this process is much less effective in ejecting them from the cluster. Accordingly, we neglect this process here.

### 2.3 Binary hardening processes

In all of the dynamical encounter processes considered in this work, viz., collisional hardening (described in this subsection), and dynamical formation and destruction processes (described in the next subsections), we shall assume the orbits to be circular, i.e., neglect their eccentricity. This is, again, a simplification used for a first look. However, it is well-known from extensive numerical simulations that a large majority of the binaries formed by tidal capture are circular or nearly so (Portegies Zwart et.al., 1997b), due to the rapid circularization which follows capture. Since our main concern here is with dynamically-formed binaries, this approximation may well be a reasonable one for describing overall evolutionary properties of such binary populations.

#### 2.3.1 Hardening in pre-X-ray binary (PXB) phase

As explained in detail in Banerjee & Ghosh (2006), referred to henceforth as BG06, the processes that harden binaries are of two types, viz., (a) those which operate in isolated binaries, and are therefore always operational, and (b) those which operate only when the binary in a globular cluster. In the former category are the processes of gravitational radiation and magnetic braking, and in the latter category is that of collisional hardening. As discussed in detail in BG06, collisional hardening, which increases with increasing , dominates at larger orbital radii, while gravitational radiation and magnetic braking, which increase steeply with decreasing , dominate at smaller orbital radii. It is these processes that harden a compact binary from its pre-X-ray binary (PXB) phase, during which its orbit is still not narrow enough for the companion (mass donor) star to come into Roche lobe contact, to the state where this Roche lobe contact does occur, at which point the companion starts transferring mass to the degenerate star, and the system turns on as an X-ray binary (XB) — either a CV or a LMXB, depending on the nature of the degenerate accretor.

Consider gravitational radiation first. The relative angular momentum loss rate due to this process is:

(7) |

Here, as before, is the mass in solar units of the degenerate primary (neutron star or white dwarf) which emits X-rays when accretion on it occurs during the mass-transfer phase of the compact binary, is the mass of its low-mass companion in solar units, and the unit of the binary orbital radius is the solar radius. We shall use these units throughout the work.

Now consider magnetic braking. The pioneering Verbunt-Zwaan (Verbunt & Zwaan, 1981) prescription for this process has been reassessed and partly revised in recent years, in view of further observational evidence on short-period binaries available now (for further details, see discussions in BG06 and references therein), and modern prescriptions are suggested in van der Sluys et.al. (2005). From these, we have chosen for this work the following one which preserves the original Verbunt-Zwaan scaling, but advocates an overall reduction in the strength of the magnetic braking process:

(8) |

Here, is the radius of the companion. Note that the strength of magnetic braking is still a matter of some controversy; while the evidence cited in the above reference argues for a reduction from the original value, it can also be argued that the presence of the well-known “period gap” in the period distribution of CVs requires a strength comparable to the original one. We have adopted here a recent prescription which is reasonably simple and adequate for our purposes: our final results do not depend significantly on the strength of this process.

Consider finally collisional hardening. As indicated earlier, it is a stochastic process, for whose continuous limit we use the prescription of Shull (1979), as has been done previously in the literature (see BG06 for a discussion). According to this prescription, the rate of increase of orbital binding energy of a compact binary due to collisional hardening is given in this limit by:

(9) |

Here, is the mass of the stars in the static background representing the cluster. We shall use and km sec as the units of and respectively. In the above units, the value of for Galactic globular clusters typically lie between and (BG06). The relation between and is:

(10) |

and the angular momentum loss rate is related to the shrinkage rate of the orbit , or hardening, as:

(11) |

The and terms on the right-hand side of Eqn. (11) are nonzero during mass transfer in the XB phase. In the PXB phase, , so that is simply related to as (see BG06 and references therein):

(12) |

Using Eqns. (10) and (9), we have in this case,

(13) |

The total rate of loss of orbital angular momentum due to the above three processes is:

(14) |

#### 2.3.2 Hardening in X-ray binary (XB) phase

As mass transfer starts upon Roche lobe contact, its effect on the angular momentum balance in the XB must be taken into account, in the manner described below. Note first that, for the radius of the Roche-lobe of the companion, we can use either the 1971 Paczyński approximation:

(15) |

which holds for , or the 1983 Eggleton approximation:

(16) |

which holds for the entire range of values of the mass ratio . Both approximations have been widely used in the literature, and they give essentially identical results for the mass ratios of interest here. We have used the Paczyński approximation here for simplicity of calculation.

At the Roche-lobe contact point, must be equal to the companion radius, the value of which is for a companion of (see above), according to the mass-radius relation for low mass stars (Ghosh, 2007). For , this translates into an orbital radius of at Roche lobe contact, using Eqn. (15). After this, the companion continues to remain in Roche-lobe contact as the binary shrinks further, and continues to transfer mass (van den Heuvel, 1991, 1992). In other words, we have

(17) |

throughout the XB phase. During this phase, the binary is already narrow enough that the collisional hardening rate is quite negligible compared to those due to gravitational radiation and magnetic braking.

Since no significant mass loss is expected from the XB in this phase, we have

(18) |

Combining Eqns. (11), (17) and (18) with a mass-radius relation for the companion of the form

(19) |

we find:

(20) |

Here, is the effective total rate of loss of angular momentum, since the collisional-hardening contributions are negligible, as explained above.

For the low-mass main sequence companions that we consider here, . However, when the mass of the companion becomes less than about 0.03, it becomes degenerate, so that (Ghosh, 2007). This results in a widening of the orbit () from this point onwards, which we do not follow here, since our study is not aimed at such systems, as explained in Sec. 4.1. This change-over point is, of course, that corresponding to the well-known period minimum of minutes in the orbital evolution of CVs and LMXBs(van den Heuvel, 1992). Henceforth, we denote the value of at the period minimum by , and we terminate the distributions of and in at a minimum value of in the figures shown in this work. Thus, in Fig. 2, we display the hardening rate against , beginning from a wide PXB phase, going into Roche lobe contact, and continuing through the mass-transfer XB phase upto the above period minimum. Note that has a very weak dependence on during the XB phase, which may have interesting consequences, as we shall see later.

### 2.4 Binary formation processes

Compact binaries with degenerate primaries and low-mass companions are formed in globular cluster (henceforth GC) cores primarily by means of two dynamical processes, namely, (i) tidal capture (tc) of a degenerate, compact star (white dwarf or neutron star) by an ordinary star, and (ii) an exchange encounter (ex1) between such a compact star and a binary of two ordinary stars, wherein the compact star replaces one of the binary members. Accordingly, the total rate of formation of compact binaries per unit binary radius, , consists of the above tc rate and ex1 rate :

(21) |

where is the orbital radius of the compact binary so formed. We now consider the rates of formation by tidal capture and by exchange.

#### 2.4.1 Tidal capture

In a close encounter between a compact star of mass and an ordinary star of mass with a distance of closest approach , tidal capture can occur if their relative speed is less than an appropriate critical speed , which we discuss below. The cross section for encounters within this distance is given by the well-known form (Spitzer, 1987):

(22) |

which gives the *differential* cross section for tidal capture around
as:

(23) |

The first terms in the right-hand sides of Eqs. (22) and (23) are the obvious geometrical cross sections and the second terms are due to gravitational focusing (also see below). It is clear that the latter terms dominate when is small, as is the case for the range of values of relevant to the problem we study here. We shall return later to the actual numerical values of of interest to us in this study.

After being tidally formed, the binary is believed to circularize very rapidly to an orbital radius , assuming conservation of angular momentum (Spitzer, 1987). Accordingly, the differential cross-section in terms of is given by:

(24) |

Here, is the critical velocity in terms of , obtained by setting in Eq. (25) below.

In a sense, the whole cross-section as expressed above may be regarded
as “geometrical”, if we look upon pure considerations of
Newtonian gravity as being geometrical. Details of the essential
astrophysics enter only when we calculate the critical speed ,
and an inversion of this relation (together with other plausible
requirements; see below) then readily gives us the range of over
which tidal capture is physically admissible. This is an interesting
topic, with literature going back to the mid-1970s and earlier, and we
summarize in this section those essential points which we need in this
work. The basic physics of tidal capture is of course that, during a close
encounter, the degenerate compact star excites non-radial oscillation
modes in the normal companion star through tidal forcing (in an
encounter between two normal stars, each excites oscillations in the
other): the energy required to excite these oscillations comes from
the kinetic energy of relative motion of the two stars, so that if
enough energy is extracted from this source by exciting these modes,
the stars become bound after the encounter. This energy condition
readily translates into one between and , giving an upper
limit on velocity for a specified as above, or, as
expressed more commonly, an upper limit on the distance of closest
approach for a specified velocity (actually, often a
*distribution* of velocities, e.g., a Maxwellian, with a specified
parameter in practical situations, as we shall see below).

The above relation between and has been calculated in the literature at various levels of detail. The pioneering estimates given in Fabian et.al. (1975) or earlier works basically employ the impulse approximation for calculating the gain in the internal energy of the tidally-perturbed star, wherein the changes in the positions of the two stars during the tidal interaction are neglected. A clear account of the procedure is given in Spitzer (1987), where the final result is evaluated for two normal stars of equal masses. Upon generalizing this procedure appropriately to the problem we study, where we have (a) unequal stellar masses and , and (b) the fact that only the normal star of mass undergoes tidally-induced oscillations, we obtain the following relation between and :

(25) |

Here, is the root-mean-square radius of the companion star, i.e., its radius of gyration which is given in the polyrtopic approximation as in terms of the companion’s radius (Spitzer, 1987).

To obtain the overall rate of tidal capture in the GC core of volume per unit interval in around , we first consider this rate around a particular value of the above relative velocity of encounter, i.e., , in terms of the above differential cross-section, remembering that the rate of encounter scales with the product of the densities and of compact stars and normal stars respectively. We then average this rate over the distribution of , obtaining the form:

(26) |

where the angular brackets indicate an average over the -distribution.

For the actual averaging, we adopt in this work a Maxwellian distribution , as has been widely done in the literature. A normalized Maxwellian is

(27) |

where is the velocity dispersion introduced earlier, for which we adopt the canonical value 10 km s in the numerical calculations (also see below).

With the aid of Eqns. (24), (25) and (27), we perform the averaging and obtain:

(28) |

The terms and above arise due to what we described respectively as the geometrical term and the gravitational focusing term in the discussion below Eq. (23). Eqns. (26) and (28) together give the total tidal capture rate as:

(29) |

where is the Verbunt parameter describing the total two-body encounter rate in the cluster core, as introduced earlier, and we have ignored compared to , which is an excellent approximation for the range of or relevant here.

We show in Fig. 3 given by Eq. (29) as a function of : this tidal capture cross-section is nearly constant for , and decreases rapidly at larger . At this point, we need to invoke additional physical arguments in order to estimate the range of values of or over which tidal capture is actually possible, and use the above cross-section only over this range for our calculations. The lower bound to the above range comes from the requirement that the two stars must form a binary and not merge into each other, and the upper bound comes from the requirement introduced earlier that enough energy of relative motion between the two stars must be absorbed by the tidally-excited oscillation modes that the stars become bound. Consider the lower bound on first. Clearly, a minimum value of this bound must be the sum of the stellar radii, which in our case leads to the bound . A more conservative bound comes from the requirement that the companion must underfill its Roche lobe after the binary has formed, i.e., , which, with the aid of Eq. (15) and , yields for the masses and we have here. The idea behind the latter requirement is apparently that if the companion overfills its Roche lobe at this point, the ensuing mass transfer is likely to lead to a merger. This seems reasonable at first, but detailed N-body simulations of recent years have suggested that this requirement may, in fact, be too restrictive. In the simulations of Portegies Zwart et.al. (1997b), which included stellar evolutionary effects according to the scheme of these authors, systems which violated the latter requirement but satisfied the former one were allowed to evolve, with the result that details of the evolution determined which systems merged and which did not. In fact, these authors found a lower limit on of approximately for tidal capture with an average companion mass very similar to ours, which is to be compared with the limits from the first requirement above, and from the latter. In view of this, we have adopted the lower bound of for our calculations here, as shown in Fig. 3.

Consider now the upper bound on . We have already given the relation between and by Eq. (25) in the impulse approximation. Remembering that for a Maxwellian, the above relation yields, for a canonical value = 10 km s as given above, an upper limit of for a polytropic index and one of for . Note that these bounds of are larger than those given for two stars of equal mass (roughly 8 for and 11 for ) in Table 6.2 of Spitzer (1987) by a factor of since scales with the mass-ratio in this manner in the impulse approximation, as can be seen readily from Eq. (25), remembering that for the companions we consider here. That should increase with increasing is qualitatively quite obvious, since, other things being equal, a higher value of the mass ratio excites tidally-forced oscillations of larger amplitude. We return below to the question of the exact scaling with this mass ratio.

As has been realized long ago, the impulse approximation is of limited validity, working best when the frequency of perturbation (i.e., tidal forcing) is not very different from those of the stellar oscillation modes that are excited by this perturbation (Fabian et.al., 1975; Spitzer, 1987). Since this is not the case for the values of estimated above, we need more accurate results, which come from detailed computations of the total energy dissipated by the above excited modes. Such numerical computations were pioneered by Press & Teukolsky (1977), and detailed results were established for various situations by several groups of authors in the mid-1980s, including Lee & Ostriker (1986) and McMillan et.al. (1987), which have been extensively used since. These results have shown that the exact upper bounds on are considerably smaller than those given by the impulse approximation, as may have been expected, since the forcing frequency falls far below those of the oscillation modes at such large separations as are given by this approximation, and the efficiency of exciting these modes drops rapidly. Some exact results are given in Table 6.2 of Spitzer (1987) from the above references, but only for the equal-mass case, where the above upper bound is 2.4 for and 3.4 for .

For our purposes here, we need to obtain the above upper bounds for our
mass ratio , which we do by doing a power-law
fit of the form to the results given
for various values of the degenerate/normal star mass-ratios in Table 3
of Lee & Ostriker (1986). This yields (note that the quantity
listed in Table 3 of Lee & Ostriker (1986) is the impact parameter defined
by these authors; scales as , as shown in their
paper). The interesting point about this scaling is that it is stronger
than that given above by the impulse approximation, which corresponds to
. Clearly, then, the impulse approximation fails to extract
the entire scaling with . The reason for this appears to be
related to nonlinear effects in exciting and dissipating tidally-induced
oscillations, but needs to be investigated further^{2}^{2}2Note that
this discrepancy is even stronger for the case where both stars are
normal, main-sequence ones, since in that case, as
can be shown readily from Table 2 in the above Lee-Ostriker reference.
An obvious line of reasoning for this would be that larger nonlinear
effects may be expected when two normal stars force tidal oscillations
in each other, but we shall not speculate on this any further here..
With the above value of , the upper bound for
our mass-ratio here is 4.1 for and 5.7 for . As the latter
value of the polytropic index is believed to give a better
representation of a low-mass main-sequence companion of the kind we
are considering here, we adopt here. With
and the value of given earlier, this translates into an
upper bound on as , which we can adopt for
these calculations.

Thus we find a range of values over which tidal capture is expected to be effective in the problem we study here. Consider now how the tidal-capture cross-section is expected to fall off at the bounds of this range. At the upper bound, the cut-off is not sharp, of course, as there is a distribution of velocities. In other words, the upper bound as given above corresponds to a suitable average (actually, root-mean-square in this case) velocity, so that at any , there will be some stars in the distribution whose velocities are sufficiently below this average that tidal capture will be possible for them. Of course, their number will decrease as increases, producing a “tail” in the tidal capture cross-section whose shape is determined by that of the velocity distribution. We have used a Maxwellian distribution here, which gives the tail seen in Fig. 3, which falls off rapidly beyond . We shall use this fall-off profile in our calculations: other profiles will not make a large difference. At the lower bound, in view of the discussion given earlier, we expect the cross-section to actually fall off gradually from about to , rather than being cut off sharply at , but we shall ignore this complication here.

We close this discussion of tidal capture with some observations on the many investigations, conclusions, and points of view that the subject has now seen for more than three decades. From the pioneering suggestion and an essentially dimensional estimate of Fabian et.al. (1975), detailed calculations of the 1980s and ’90s have reached interesting, and sometimes contradictory, conclusions. For example, concerns that energy dissipation by tidally-induced modes may lead to a large distention of the companion and so to a merger have been confronted with results from detailed computations of the nonlinear damping of the primary modes by coupling to other, high-degree modes, which suggested that the damping took place far more rapidly than thought before, and the energy dissipated was too small to have a significant effect on the companion’s structure. We here have a adopted a somewhat moderate view that tidal capture is plausible, but efficient over only a restricted range of or . This view is supported by (a) recent observational demonstration that the number of X-ray sources in Galactic globular clusters scale with their Verbunt parameter , i.e., the two-body encounter rate (Pooley et al., 2003), as described earlier, and (b) recent N-body simulations of Portegies Zwart et.al. (1997b) showing tidal capture over a considerable range of , admittedly under the algorithms adopted by these authors. Consider, finally, our suggested range of radii for efficient tidal capture, , as given above, in the context of other suggested ranges. Values in the range have been thought plausible by Podsiadlowski et.al. (2002), while Portegies Zwart et.al. (1997b) have demonstrated tidal capture over a range . We here advocate a range (depending on ), which is between the two, and still quite modest.

#### 2.4.2 Formation by exchange

Exchange encounters between binaries and single stars with arbitrary mass ratios has been extensively studied by Heggie, Hut & McMillan (1996). They performed detailed numerical scattering experiments, using the automatic scattering tools of the STARLAB package. From the resulting exchange cross sections, they obtained a semi-analytic fit of the form:

(30) |

Here, is the orbital radius of the initial binary, is the mass of the escaping star, is the companion mass, is the mass of the incoming star, and . is the dimensionless cross section which is a function of these masses only and which is given by Eq. (17) of Heggie, Hut & McMillan (1996). We use Eqn. (30) to obtain the cross sections for the exchange process ‘ex1’ described above, but one essential point needs to be clarified first.

The radius of the compact binary formed by exchange is not the
same as the radius of the original binary undergoing exchange.
Therefore, a relation between and is required, since in
Eqn. (30) represents the radius of the initial
binary, *not* the radius of the compact binary formed by
exchange. According to the binary-hardening rule of Heggie (Heggie, 1975),
the final compact binary must, on an average, be harder, i.e., have a
larger binding energy. We performed illustrative scattering experiments
with circular binaries and incoming stars with mass ratios of interest
to us in this study, using the scattering tools of STARLAB. The resulting
distribution of the change in orbital radius is shown in
Fig. 4, and is seen to be highly asymmetric.

The long tail towards implies that the binary radius
*increases* in many scatterings. This does not of course contradict
the above Heggie rule, since the increase of mass due to exchange (the
mass of the incoming compact star, 1.4, is a factor
times the mass of the outgoing low-mass star, 0.6) increases the
binding energy by itself by the above factor. From these experiments, we
see that the peak of the distribution corresponds to a shrinkage of the
binary by about 25 per cent. On the other hand, the *average*
change in binary radius, calculated from the above distribution, is much
closer to zero due to the above long tail of the distribution on the
side, so that we can take for our
purposes here without much error.

The total Maxwellian-averaged rate of formation of compact binary by this type of exchange (ex1) in the GC core is then:

(31) |

Here, is the distribution function of the orbital radii of the primordial stellar binaries in the cluster core. For primordial binaries, we can take the widely-used distribution (i.e., a uniform distribution in ) (Kraicheva et.al., 1978), with a lower bound at , corresponding to the smallest possible radius for a binary of two main-sequence stars. The ex1 rate is shown in Fig. 3.

### 2.5 Binary destruction processes

A compact binary can be destroyed by two major processes. First, an
encounter with a star which has a relative speed higher than an
appropriate critical speed (Hut & Bahcall, 1983) can lead to its dissociation
(dss). Second, in an exchange encounter (ex2) of this binary with a
compact star, the latter can replace the low-mass companion in the
binary, forming a double compact-star binary consisting of two neutron
stars, two white dwarfs, or a neutron star and a white dwarf, all with
masses . This, in effect, destroys the binary as an
X-ray source (as accretion is not possible in such a system), and so
takes it out of our reckoning in this study. This is so because such a
system is not an X-ray source, and it is essentially impossible for one
of the compact stars in such a system to be exchanged with an ordinary
star in a subsequent exchange encounter, since is much
lighter than . The total destruction rate D(a) *per
binary* is thus the sum of the above dissociation and exchange rates:

(32) |

We now discuss the rates of these two processes.

#### 2.5.1 Dissociation

To estimate the dissociation rate of compact binaries, we use the
results of scattering experiments of Hut & Bahcall (1983). The
Maxwellian-averaged dissociation rate (dss) *per compact binary*
is then given by

(33) |

From Hut & Bahcall (1983), we adopt

(34) |

a relation which was obtained by these authors by fitting the results of their scattering experiments with analytical models. Here, is the threshold relative velocity for ionization (see Sec. 2.5), given by:

(35) |

As these authors pointed out, Eqn. (34) is an asymptotic form, which works well only for significantly hard binaries, i.e., those with . This condition is of course satisfied for the compact binaries that we are interested in here.

We show in Fig. 3 the above dissociation rate, whose essential variation with is seen by combining Eqs. (34) and (35), which yields the form , where is a constant. Thus, the dissociation rate is quite negligible for , reflecting the fact that it is essentially impossible to dissociate very hard binaries. As increases, the rate rises extremely sharply at first (the initial rise is determined by the exponential), and eventually scales as for .

#### 2.5.2 Destruction by exchange

By arguments similar to those given in Sec. 2.4.2, we arrive at a
Maxwellian-averaged rate of this type of exchange (ex2) *per
compact binary* which is:

(36) |

and which is also shown in Fig. 3. This rate scales with simply as . Note the different magnifications used for different curves in Fig. 3 in order to make all of them clearly visible. Of the two destruction processes, dominates completely at all orbital radii of interest in our study (reflecting the fact that dynamically-formed binaries in GC cores are so hard that they cannot be dissociated or “ionized” by further encounters in that GC core), but the fast-rising eventually overtakes it at , corresponding to very soft binaries.

### 2.6 The numerical method

Equation (5) for the evolution of compact binary populations is a partial differential equation (PDE) of hyperbolic type, with similarities to wave equations, as pointed out earlier. We solved this equation using a Lax-Wendorff scheme (Press et.al., 1992). This involves dividing the range of and in a discrete mesh () of constant space intervals () and time intervals (). The PDE is then discretised into a set of linear difference equations over this mesh, which is solved numerically.

We denote by the value of at the th time step and the th point in . Discretisation of Eqn. (5) according to the Lax-Wendorff scheme is a two-step process:

(37) |

For a chosen mesh-interval , Eqn. (37) will be numerically
stable only if the time-step is chosen to be small enough that it
obeys the *Courant condition* (Press et.al., 1992) throughout the mesh:

(38) |

where is the maximum value of within the -range of the mesh.

We chose Lax-Wendorff scheme among the various existing schemes for solving hyperbolic PDEs primarily because it appears to be the only explicit method that does not have any significant numerical dissipation (Press et.al. 1992, and references therein) and is at the same time numerically stable, provided that the time step is chosen according to the Courant condition. This point is important, since numerical dissipation can significantly affect the computed evolution of and the X-ray binary population, as we observed while trying other methods, e.g., the so-called staggered-leapfrog method. Other instabilities, e.g., the mesh-drifting instability (Press et.al., 1992), also appeared to be insignificant in the method we chose.

## 3 Results

### 3.1 Evolution of compact-binary distribution

A typical result from our computed evolution of the compact-binary distribution function is shown in Fig. 5, wherein the surface is explicitly displayed in three dimensions. The GC parameters chosen for this run were , and , similar to those of the well-known Galactic cluster 47 Tuc. The distribution function is seen to evolve as a smooth surface, with the compact binary population growing predominantly at shorter radii (, say). We start with a small number of binaries at following various distributions, and find that the distribution at large times Gyr is quite independent of these initial conditions, being determined entirely by the dynamical processes of formation and destruction, and by the various hardening processes detailed earlier. Note that, since the point of Roche lobe contact corresponds to in our study, as explained earlier, that part of the distribution which is shortward of this radius corresponds to XBs, while that part longward of it corresponds to PXBs.

To further clarify the nature of this evolution, slices through the above surface at various points along time axis and -axis are shown in Figs. 6 and 7, in the former figure the abscissa being also marked in terms of the orbital period , readily calculable in terms of and the stellar masses with the aid of Kepler’s third law, assuming conservative mass transfer during the XB phase. Figure 6 shows that increases with time, roughly preserving its profile for Gyr or so. This profile consists of a roughly uniform distribution in for short orbital radii, , say, corresponding to roughly, and a sharp fall-off at larger radii and orbital periods. Figure 7 shows that at a given increases with time and approaches saturation on a timescale Gyr or so, this timescale being longer at at smaller values of .

Figures 6 and 7 suggest that a regime of roughly self-similar evolution may be occurring in our model binary population at times beyond 1 Gyr or so, in the following way. An asymptotic profile of is established on the timescale of a 1 Gyr or so, which thereafter evolves roughly self-similarly towards a saturation strength on a timescale Gyr or so. We shall discuss the origins of such behavior in detail elsewhere, since, as explained in Sec. 4.1, our model of orbital evolution requires additional ingredients before it can be compared with observations of X-ray binaries. However, the following qualitative remarks are appropriate here.

First, the origins of the establishment of the above self-similar profile in a Gyr or so (independent of the initial distribution we start from) clearly lie in the two terms that describe binary formation and hardening on the right-hand side of Eq. (5), namely, and respectively. The latter term can be written qualitatively in the form , where is the overall hardening timescale, which is well-known to be of the order of a Gyr or so (see BG06 and references therein). This timescale, which is also that on which a given binary passes from the large- end of the distribution shown in Fig. 5 to the small- end, is obviously the timescale that establishes the above profile. The shape of this profile, as detailed above, seems related to those of the tidal-capture rate (see Fig. 3) and the hardening rate (see Fig. 2). In particular, note that the former rate is roughly constant over , and the latter roughly so for .

Second, the subsequent, roughly self-similar evolution of the above profile occurs on a (longer) timescale whose origins clearly lie in the binary destruction term on the right-hand side of Eq. (5), namely, , since this term can be cast in the qualitative form , where is the saturation time Gyr. Whereas the earlier term describes the passage or “current” of binaries through the distribution, as described earlier, the term becomes important as increases, preventing from becoming arbitrarily large by enforcing saturation at the point where the rates of formation and destruction balance. As scales with , as shown above, and , we expect saturation to occur at earlier times at larger radii, as seen in Fig. 7.

### 3.2 Number of X-ray binaries in globular clusters

The total number of X-ray binaries in a GC at any time can be computed directly from our approach by integrating over the range of relevant for XBs, viz., , where is the value of corresponding to the period minimum minutes, and is the value of at the first Roche lobe contact and onset of mass transfer, as explained earlier. We have:

(39) |

Taking an evolutionary time Gyr as representative, we can therefore determine at this point in time, and study its dependence on the Verbunt parameters and that describe the essential dynamical properties of globular clusters in this context, as explained earlier. By doing so, we can attempt to make qualitative contact with the systematics of those recent observations of X-ray binaries in globular clusters which we have described earlier (Pooley et al., 2003). To this end, we computed values of over a rectangular grid spanning over and . (Of course, not all the points on the grid would be directly relevant for comparison with observation, since the observed globular clusters lie only along a diagonal patch on this grid, as shown in Fig. 1. However, in this introductory study, we wished to establish the theoretically expected trends of variation with and , and so performed computations of over the entire rectangular grid)

For a specified grid point, i.e., a pair of values of the Verbunt parameters, we obtained values of , and with the aid of Eqs. (1), (2) and the virialization condition:

(40) |

which were used for the computation at this grid point. We chose this
prescription for the sake of definiteness, because values of
are not known, in general, at a computational grid point, without
which a pair of Verbunt parameters cannot specify all three variables
, and . This also introduced a certain uniformity
of treatment of all grid points, which, we thought, would clarify the
theoretically expected trends. On the other hand, this did lead to a
feature at high values of and low values of ,i.e., in
that part of the grid which is completely devoid of observed globular
clusters at this time (and which, in fact, may *actually* contain
no clusters, because such combinations of and may
not be possible in nature), which appears unphysical, as we discuss
below. Observationally, we know, of course, that some clusters appear
fairly virialized and some do not, but any spread in applied
over the grid points would have been arbitrary, and would have led to
a scatter, masking the systematic theoretical trends without purpose.
Finally, throughout these computations, we used representative values
for (a) the primordial binary fraction , namely, 10 percent, and
(b) compact star fraction , namely, 5 percent.

Figure 8 shows the computed surface . There appears to be a “fold” in this surface, in a direction roughly parallel to the axis, located around . From this fold, if we go towards higher values of , then, for any given value of , decreases with increasing . This is a signature of the compact-binary destruction processes detailed in the previous section, whose strengths increase with increasing . Thus, the above value of corresponding to the fold seems to be a good estimate of the threshold above which these destruction processes dominate. At constant , the variation with is quite straightforward: simply increases monotonically with increasing , reflecting the fact that the formation rates of compact binaries, as described in the previous section, increase with increasing .

To further clarify these trends, and to facilitate comparison with those obtained from the “toy” model of BG06, we display in Fig. 9 vs. , as was done in that reference. The motivation is as follows. It was shown in BG06 that the toy model of these authors leads to the scaling that was a function of alone, which was a monotonically increasing function of , for which the toy model gave a very simple, analytic form. Our purpose in Fig. 9 is to see how much of this scaling survives the scrutiny of a more detailed model, such as presented here. As the figure shows, this scaling does carry over approximately, although some details are different. is still almost a function of alone (except at the very highest values of ), showing that this scaling of the toy model carries over approximately to more detailed ones, thereby giving an indication of the basic ways in which dynamical binary formation and destruction processes work. The above “universal” function of is, except for a feature at low values of which we discuss below, still a monotonically increasing one, reflecting the increasing strength of dynamical binary-destruction processes with increasing . However, the shape of the function is different in detail now, as may have been expected.

We now discuss the low- feature referred to above: at the lowest values of , seems to rise again, reflecting an apparent drop in . This is difficult to understand, since binary-destruction effects are negligible at these values of . Actually, this is an artifact of the way in which we fixed the essential cluster parameters , and from specified values of the Verbunt parameters for the computational grid (as explained above), which can be seen as follows. With the assumption of virialization, as done for this purpose, the velocity dispersion can be expressed in terms of the Verbunt parameters in a manner analogous to that used in Eq. (3), the result being . This relates to , so that the latter influences the Maxwellian-averaging process involved in the calculation of the tidal capture cross-section described in Sec. 2.4.1, since the parameter of the Maxwellian then scales as . At small values of , becomes small, which reduces the tidal-capture rate, as Eq. (29) readily shows. This is completely unphysical, of course, since has nothing to do physically with the tidal capture rate. Rather, it is an artifact produced by the way we (artificially) related to for computational convenience. Accordingly, we ignore this low- feature in all further considerations.

## 4 Comparison with Observation

### 4.1 Applicability of our study

Before attempting to compare our results with observations, we review in brief some essential ingredients of our model study at this stage, so as to clarify which of our results can be so compared, and which need inclusion of further components before this can be meaningfully done. A major ingredient that is incomplete at this stage is our description of the orbital evolution of the binary, since it neglects nuclear evolution of the low-mass companion star altogether. While this may not be very unreasonable for CV systems or for short-period LMXBs with orbital periods between hours and the above period minimum of minutes, it is completely inadequate for other LMXBs, where the stellar evolution of the companion plays a crucial role, which has been studied by many authors. In particular, recent studies by Podsiadlowski et.al. (2002) and Pfahl et.al. (2003) have demonstrated the large range of possibilities covered by such evolution with realistic stellar evolutionary codes, performing a Monte Carlo binary population synthesis study in the latter reference with the aid of the library of evolutionary sequences described in the former. We plan to include stellar evolutionary effects in a subsequent work of the series and are assessing various methods of doing so. One possibility is to start with a semi-analytic scheme along the lines of the SeBa model as described in Portegies Zwart et.al. (1997b), and to continue with a semi-analytic approximation to a more elaborate library of evolutionary sequences, such as described above.

Since most of the XBs in the Galactic GC data of Pooley et al. (2003) are
CVs, our scheme should be able to describe the overall properties of
these XB populations reasonably well. Even so, we shall make no
attempt here to compare our results on orbital period distribution
with the observed CV distribution, since the CVs in the latter
distribution are almost exclusively from outside globular clusters,
where dynamical formation is not relevant. We have here recorded the
orbital-period distribution that comes from our computations (at
this stage) only as a natural intermediate step. It can perhaps be
compared with observation when the orbital-period distribution of
CVs *in GCs* becomes observationally established. For LMXBs,
where the observed orbital-period distribution at this time also
consists overwhelmingly of those outside GCs, there is of course no
question of comparison at this stage, for the reasons given above.
Thus, our main aim here is to put in the observational context our
results on the numerical properties of XB populations in GCs in
relation to the GC parameters.

### 4.2 Ultracompact X-ray binaries

In recent years, a subset of LMXBs in GCs, in the Milky Way and possibly also in elliptical galaxies, have received much attention because of (a) their high, persistent brightness ( erg s), which would make them dominate the high end of the luminosity functions of X-ray binaries in ellipticals (Bildsten & Deloye, 2004) and (b) their very close orbits with hr or so, sometimes as short as minutes, the classic example being the 11 min binary 4U 1820-30 the Galactic cluster NGC 6624. These are the ultracompact X-ray binaries (henceforth UCXBs), which are thought to consist of neutron stars in ultracompact orbits with very low-mass degenerate dwarf companions () as mass donors. The evolutionary origin of UCXBs is of much current interest, and proposals for such origins include (a) direct collisions between red giants and neutron stars in GC cores, as a consequence of which the red-giant envelope can either be promptly disrupted (Ivanova et.al., 2005) or be expelled more slowly in a common-envelope phase, and (b) usual LMXB evolution with the initial orbital period below the “bifurcation period” of about 18 hrs (Podsiadlowski et.al., 2002). A natural point that arises, therefore, is about the role of UCXBs in our study, and the general importance of the above channels of formation in relation to the ones we have described above, which we now consider in brief.

The key feature of UCXBs from the point of view of our study is that the number of UCXBs is a tiny fraction of the total number of XBs in a GC, and so of little importance as far as is concerned. This is a general, robust feature, which follows from the basic point that the UCXBs are extremely short-lived because of their extreme brightness, so that is small at any given epoch despite their considerable birthrate. To see this in more detail, consider the UCXB birthrate of about one every year per of the mass of a GC, as given by (Bildsten & Deloye, 2004), which, together with their estimated lifetimes of years, yields an estimate of in a GC at any given time. Actually, the observed GCs in our galaxy have lower masses, in the range (Ivanova et.al., 2005). Thus a Galactic GC of like 47 Tuc will have , remembering that the birthrate scales down appropriately with the GC mass, but the lifetime remains the same. This is to be compared with the observed number of XBs in 47 Tuc of 45 (Pooley et al., 2003), which yields a fraction . We can double-check this and put it on a systematic basis with the aid of Table 1 of Ivanova et.al. (2005), wherein these authors have listed the minimum expected number of UCXBs in a number of Galactic GCs, by combining this with the total number of observed XBs obtained from Pooley et al. (2003) and other sources. For 47 Tuc, with 0.23 UCXBs and 45 XBs, the ratio is , very similar to above, and those for other sources are also similar. For example, Terzan 5 has a ratio , and NGC 6652 has a ratio .

We see from the above that UCXBs constitute such a tiny fraction of the total XB populations of Galactic GCs in terms of numbers that their effect is negligible for this work. However, in a study of the X-ray luminosity functions of GCs, their effect is expected to be crucial: if a GC contains even one UCXB, its luminosity may dominate over the combined output of all other XBs. It is the extension of this idea which has been used in recent years to argue that the luminosity function of XBs in ellipticals may be dominated by UCXBs in their GCs (Bildsten & Deloye, 2004).

### 4.3 X-ray source numbers in globular clusters

The filled squares in Fig. 8 represent globular clusters with significant numbers of X-ray binaries in them. These points generally lie near the surface in this three-dimensional plot, mostly in the vicinity of the fold described above. This is more clearly seen in the two-dimensional plot of Fig. 9, where the bulk of the observational points are indeed seen to be near the upward “knee” of the computed curves. To facilitate comparison with observations further, we show in Fig. 10 contours of constant in the (Verbunt parameters) plane. Overplotted on these are the above observed clusters (filled sqaures), where the number in the parentheses next to each indicates the total number of X-ray binaries observed in it (Pooley et al., 2003). The contours are seen to be qualitatively rather similar in shape to the curves in Fig. 9. The trend in the observed values generally follows the contours, with one exception. This is most encouraging (also see BG06) for the construction of more detailed models, and indeed rather remarkable in view of the fact that no particular attempt has been made to fit the data at this stage.

## 5 Discussion

In this paper, we have explored the results of a Boltzmann study of the evolution of compact-binary populations in globular clusters in the continuous limit, and made preliminary contacts with observations of X-ray binaries in Galactic globular clusters. Our Boltzmann approach has built into it the rates of the essential dynamical processes that occur due to star-star and star-binary encounters in dense clusters, viz., collisional hardening, binary formation by tidal capture and exchange, and binary destruction by dissociation and other mechanisms, as obtained by previous numerical studies of large numbers of such individual encounters. We stress that our Boltzmann scheme is not a Fokker-Planck one, wherein the cumulative effects of a large number of small changes in distant encounters is described as a slow diffusion in phase space. We can and do handle both small and large changes in the framework of the original Boltzmann visualization of motion through phase space (at a computational cost which is quite trivial compared to that required for N-body simulations). Indeed, the continuous limit of collisional hardening used in this paper may be looked upon as an example of a slow diffusion in -space, while some of the formation and destruction processes are examples of faster and more radical changes. Of course, all these processes are episodic in nature, and we are studying their continuous, probabilistic limit in this introductory paper. As already pointed out, Paper II will describe an explicit treatment of the stochasticity of these processes within the framework of stochastic PDEs, which the Boltzmann equation becomes under such circumstances.

### 5.1 Conclusions

We find the indications from this preliminary study to be sufficiently encouraging to attempt several steps of improvement, most of which we have already indicated in the previous sections. To recapitulate briefly, we need to provide an appropriate description of the stochastic processes, which we do in Paper II. We need to introduce a mass function for the background stars in the globular cluster core, and handle non-circular orbits formed in the encounter processes. We need to assess the possible importance of binary-binary interactions in this problem, which we have neglected altogether so far. We need to include essential aspects of stellar evolution of the companion in our orbital-evolution scheme, particularly for LMXBs. In a more ambitious vein, we need to consider the evolution of the stellar background representing the cluster core, which we do in Paper III. As the core collapses, the collapse stalls due to binary heating, and possible gravothermal oscillations occur, the core parameters and evolve appropriately, and so do the Verbunt parameters and . Binary-population evolution with such evolving GC parameters is an interesting problem in itself, even if we do not explicitly consider the back reaction of binary evolution on the evolution of its background.

The scaling of with the two Verbunt parameters we already found here seems to be among the basic building blocks of our understanding of how globular clusters cook up their gross overabundance of X-ray binaries through an interplay between dynamical formation and destruction. It remains to be seen if there are other such building blocks which have not been investigated so far.

## References

- Banerjee & Ghosh (2006) Banerjee S. and Ghosh P., 2006, MNRAS, 373, 1188. (BG06)
- Bildsten & Deloye (2004) Bildsten L. and Deloye C. J., 2004, ApJ, 607, L19.
- Fabian et.al. (1975) Fabian A. C., Pringle, J. E. and Rees, M. J. 1975, MNRAS, 172, 15P.
- Gao et.al. (1991) Gao B., Goodman J., Cohn H. and Murphy B., 1991, ApJ, 370, 567.
- Ghosh (2007) Ghosh P., 2007, Rotation and Accretion Powered Pulsars, World Scientific Publications.
- Harris (1996, revised in 1999) Harris W.E., 1996, Globular Clusters in the Milky Way, VizieR On-line Data Catalog: VII/195.
- Heggie (1975) Heggie D., 1975, MNRAS, 173, 729.
- Heggie, Hut & McMillan (1996) Heggie D., Hut P. and McMillan S.L.W., 1996, ApJ, 467, 359.
- Heggie & Hut (2003) Heggie D. and Hut P., 2003, The Gravitational Millon-Body Problem, Cambridge University Press.
- Hut (2001) Hut P., 2001, in Astrophysical Supercomputing Using Particles, IAU Symposium, Vol. 208, eds. Makino J. and Hut P.
- Hut & Bahcall (1983) Hut P. and Bahcall J.N., 1983, ApJ, 268, 319.
- Hut, McMillan & Romani (1992) Hut P., McMillan S. and Romani R.W., 1992, ApJ, 389, 527.
- Hut et.al. (1992) Hut P. et. al., 1992, PASP, 104, 981.
- Ivanova et.al. (2005) Ivanova N. et. al., 2005, ApJ, 621, L109.
- Kraicheva et.al. (1978) Kraicheva Z.T., Popova E.I., Tutukov A.V. and Yungelson L.R., 1978, in “Nonstationary Evolution of Close Binaries”, ed. Zytkow A., Polish Academy of Sciences, p.25
- Lee & Ostriker (1986) Lee H.M. and Ostriker J.P., 1986, ApJ, 310, 176.
- Liu et.al. (2001) Liu Q.Z., van Paradijs J. and van den Heuvel E.P.J., 2001, A&A, 368, 1021.
- Makino (1996) Makino J., 1996, ApJ, 471, 796.
- Makino & Aarseth (1992) Makino J. and Aarseth S.J., 1992, PASJ, 44, 141.
- Makino & Taiji (1998) Makino J. and Taiji M., 1998, Scientific Simulations with Special-Purpose Computers–the GRAPE Systems, Wiley-VCH, New York.
- McMillan et.al. (1987) McMillan S.L.W., McDermott P.N. and Taam R.E., 1987, ApJ, 318, 261.
- Pfahl et.al. (2003) Pfahl E.D., Rappaport S. and Podsiadlowski P. 2003, ApJ, 597, 1036.
- Podsiadlowski et.al. (2002) Podsiadlowski P., Rappaport S. and Pfahl E.D. 2002, ApJ, 565, 1107.
- Pooley et al. (2003) Pooley D. et.al., 2003, ApJ, 591, L131.
- Portegies Zwart et.al. (1997a) Portegies Zwart S. F., Hut, P. and Verbunt F., 1997a, A&A, 328, 130.
- Portegies Zwart et.al. (1997b) Portegies Zwart S.F., Hut P., McMillan S.L.W. and Verbunt F., 1997b, A&A, 328, 143.
- Press & Teukolsky (1977) Press W.H., and Teukolsky S.A., 1977, ApJ, 213, 183.
- Press et.al. (1992) Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., 1992, Numerical Recipes in C, Cambridge University Press.
- Shull (1979) Shull J.M., 1979, ApJ, 231, 534.
- Sigurdsson & Phinney (1993) Sigurdsson S. and Phinney E.S., 1993, ApJ, 415, 631.
- Sigurdsson & Phinney (1995) Sigurdsson S. and Phinney E.S., 1995, ApJS, 99, 609.
- Spitzer (1987) Spitzer L. Jr., 1987, Dynamical Evolution of Globular Clusters, Princeton Univ. Press.
- Sugimoto & Bettwieser (1983) Sugimoto D. and Bettwieser E., 1983, MNRAS, 204, 19p.
- van den Heuvel (1991) van den Heuvel E.P.J., 1991, in “Neutron Stars: Theory and Observation”, eds. Ventura J. and Pines D., Kluwer, Dordrecht, p.171.
- van den Heuvel (1992) van den Heuvel E.P.J., 1992, in “X-Ray Binaries and Recycled Pulsars”, eds. van den Heuvel E.P.J. and Rappaport S.A., Kluwer, Dordrecht, p.233.
- van der Sluys et.al. (2005) van der Sluys M. V., Verbunt F. and Pols, O.R. 2005, A&A, 440, 973.
- Verbunt (2002) Verbunt F., 2002, New horizons in globular cluster astronomy, ASP conf. series., 296, 245, eds. G. Piotto et al., Astron. Soc. Pacific, San Francisco.
- Verbunt (2006) Verbunt F., 2006, Highlights of Astronomy, Volume 14, Proc. XXVIth IAU General Assembly, Prague 2006, ed. van der Hucht K.A., IAU Publ, Paris.
- Verbunt & Zwaan (1981) Verbunt F. and Zwaan C., 1981, A&A, 100, L7.