text:
# Grand canonical ensemble simulation studies of polydisperse fluids

## I Introduction

Statistical mechanics was originally formulated to describe the properties of systems of identical particles such as atoms or small molecules. However, many materials of industrial and commercial importance do not fit neatly into this framework. For example, the particles in a colloidal suspension are never precisely identical to one another, but have a range of radii (and possibly surface charges, shapes etc). This dependence of the particle properties on one or more continuous parameters is known as polydispersity.

To process a polydisperse colloidal material, one needs to know its phase behavior, i.e. the conditions of temperature and pressure under which a given structure is thermodynamically stable. The main obstacles to gaining this information arise from the effectively infinite number of particle species present in a polydisperse system. Labeling these by the continuous polydispersity attribute $`\sigma `$, the state of the system must be described by a density distribution $`\rho (\sigma )`$, rather than a finite number of density variables. The phase diagram is therefore infinite dimensional, a feature that poses serious problems to experiment and theory alike.

The chief difficulty faced in experimental studies of polydisperse systems is that the infinite dimensionality of the phase diagram precludes a complete mapping of the phase behavior. Instead one is forced to focus attention on particular low dimensional manifolds (slices) of the full diagram. Typically this involves determining the system properties along some desired trajectory through the space of $`\rho (\sigma )`$. Such a strategy is often pursued in experiments on colloidal suspensions , where the phase behavior is studied along a so-called dilution line. The experimental procedure for tracking this line involves adding a prescribed quantity of colloid of some known degree of polydispersity to a vessel of fixed volume $`V`$, the remaining volume being occupied by a solvent. The number distribution of colloidal particles $`N(\sigma )`$ determines the density distribution of the suspension, $`\rho (\sigma )N(\sigma )/V`$. Since in a given substance, the relative proportions of the number of particles of each $`\sigma `$ are fixed, changing the amount of colloid added simply alters the scale of $`\rho (\sigma )`$, not its shape. Thus, by varying $`N(\sigma )`$ at fixed $`V`$ (or vice versa) one traces out a locus in the phase diagram in which only the overall scale of $`\rho (\sigma )`$ changes.

As regards theoretical studies of phase behavior, these typically endeavor to calculate the system free energy as a function of a set of density variables. The difficulty in achieving this for a polydisperse system is that the free energy $`f[\rho (\sigma )]`$ is a functional of $`\rho (\sigma )`$, and therefore itself occupies an infinite dimensional space. This renders intractable the task of identifying phase boundaries. Recently, however, progress has stemmed from the observation that it is possible to approximate the full free energy by a so-called “moment free energy” containing the full ideal gas contribution plus an excess part that depends only on certain principal moments of the full excess free energy . Doing so reduces the problem to a finite number of density variables and allows calculation of phase coexistence properties within a systematically refinable approximation scheme. Additionally the theory delivers (for the given free energy) exact results for the location of spinodals, critical points and the cloud and shadow curves. Use of this approach promises to enhance significantly our understanding of phase behavior in polydisperse systems.

In view of the approximations inherent in theoretical approaches, it is natural to consider deploying computer simulation to study the phase behavior of polydisperse colloidal fluids. Although simulations (like experiment) are restricted to studying limited regions of the phase diagram such as the dilution line, they have the advantage that can be used to investigate the same model systems as studied theoretically. Furthermore, they deliver (modulo finite-size effects) essentially exact results, providing invaluable benchmarks against which to test theoretical predictions. Sometimes too, the physical insight gleaned from simulations serves as the impetus for fresh theoretical advances.

One simulation approach for obtaining the thermodynamic properties of a polydisperse system is to simply mimic the experimental procedure. This can be achieved by employing a canonical ensemble (CE) formalism wherein a simulation box of fixed volume $`V`$ is populated by a prescribed number of particles $`N`$, whose sizes are distributed according to the desired $`N(\sigma )`$. It practice, however, it transpires that the CE represents a far from optimal framework for simulating polydisperse fluids. The principal difficulty lies with the limited range of computationally accessible particle numbers which, in any simulation, is typically many orders of magnitude smaller than found in an experiment. The resulting finite-size effects are particularly pronounced in a CE simulation because the specific realization of the disorder $`N(\sigma )`$ is fixed. This suppresses large scale fluctuations in $`\rho (\sigma )`$ and potentially leads to sampling deficiencies . Additionally, the CE suffers other drawbacks familiar from simulation studies of monodisperse fluids. For instance, relaxation times are extended because density fluctuations decay solely via diffusion; there is no direct access to information on chemical potentials; metastability and hysteresis hinder the study of phase transitions.

Experience with the simulation of monodisperse fluids has shown that use of the grand canonical ensemble (GCE) is highly effective in circumventing many of the aforementioned problems associated with the CE . As we shall show, its application in the context of polydisperse fluids retains many of the benefits of the monodisperse case. Moreover, it provides the key to improved sampling of the density distribution $`\rho (\sigma )`$. This is because within the GCE framework $`N(\sigma )`$ fluctuates as a whole, thereby capturing fluctuations in $`\rho (\sigma )`$ on all simulation length scales. Notwithstanding these advantages, however, the GCE might appear, at first sight, unsuitable for the purpose of traversing a particular trajectory through the space of the density distribution $`\rho (\sigma )`$. This is because $`\rho (\sigma )`$ ostensibly lies out-with the direct control of the simulator, its form being instead determined by the imposed chemical potential distribution $`\mu (\sigma )`$. Nevertheless, it turns out to be possible to tune $`\mu (\sigma )`$ within a histogram extrapolation scheme, in such a way as to realize a specific desired form of $`\rho (\sigma )`$. We shall demonstrate that this dual use of the GCE and histogram extrapolation methods permits a chosen phase space path to be followed efficiently.

The layout of our paper is as follows. In sec. II A we formulate the statistical mechanics for a polydisperse fluid within the grand canonical ensemble. We then describe (sec. II B) the combined GCE plus histogram extrapolation methodology for tracking a particular path through the space of $`\rho (\sigma )`$. In section III we apply the method to the problem of obtaining the dilution line properties of three types of size-disperse hard sphere fluids. The chemical potential distribution of these fluids is determined as a function of volume fraction and the results compared with the predictions of two commonly used equations of state. Finally in section IV, we discuss our findings and their implications.

## II Method

### A Statistical Mechanics

We consider a classical fluid of polydisperse particles confined to a volume $`V=L^d`$. The system is assumed to be thermodynamically open, so that the particle-number distribution $`N(\sigma )`$ is a statistical quantity. The associated grand canonical partition function takes the form:

$$𝒵_V=\underset{N=0}{\overset{\mathrm{}}{}}\frac{1}{N!}\underset{i=1}{\overset{N}{}}\left\{_V𝑑\stackrel{}{r}_i_0^{\mathrm{}}𝑑\sigma _i\right\}\mathrm{exp}\left(\beta _N\left(\{\stackrel{}{r},\sigma \}\right)\right)$$
(1)
with

$$_N\left(\{\stackrel{}{r},\sigma \}\right)=\mathrm{\Phi }\left(\{\stackrel{}{r},\sigma \}\right)\underset{i=1}{\overset{N}{}}\mu (\sigma _i).$$
(2)
Here $`N`$ is the overall particle number, while $`\beta =(k_BT)^1`$ and $`\mu (\sigma )`$ are respectively the prescribed inverse temperature and chemical potential distribution. $`\{\stackrel{}{r},\sigma \}`$ denotes the configuration, i.e. the complete set $`(\stackrel{}{r_1},\sigma _1),(\stackrel{}{r_2},\sigma _2)\mathrm{}(\stackrel{}{r_N},\sigma _N)`$ of particle position vectors and polydisperse attributes. The corresponding configurational energy $`\mathrm{\Phi }\left(\{\stackrel{}{r},\sigma \}\right)`$ is assumed to reside in a sum of pairwise interactions

$$\mathrm{\Phi }\left(\{\stackrel{}{r},\sigma \}\right)=\underset{i<j=1}{\overset{N}{}}\varphi (\stackrel{}{r}_i\stackrel{}{r}_j,\sigma _i,\sigma _j),$$
(3)
where $`\varphi `$ is the pair potential.

The instantaneous number distribution is defined by

$$N(\sigma )\underset{i=1}{\overset{N}{}}\delta (\sigma \sigma _i),$$
(4)
with $`\sigma `$ the continuous polydispersity attribute and $`N=N(\sigma )𝑑\sigma `$. We shall be concerned with the fluctuations in the associated density distribution

$$\rho (\sigma )N(\sigma )/V,$$
(5)
and the configurational energy density

$`u\mathrm{\Phi }(\{\stackrel{}{r},\sigma \})/V.`$The statistical behavior of these observables is completely described by their joint probability distribution

$$p_V[\rho (\sigma ),u]=\delta \left(uV^1\mathrm{\Phi }(\{\stackrel{}{r},\sigma \})\right)\underset{\sigma }{}\delta \left(\rho (\sigma )V^1N(\sigma )\right),$$
(6)
or more explicitly

$`p_V[\rho (\sigma ),u]`$ $`=`$ $`{\displaystyle \frac{1}{𝒵_V}}{\displaystyle \underset{N=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{N!}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle _V}𝑑\stackrel{}{r}_i{\displaystyle _0^{\mathrm{}}}𝑑\sigma _i\right\}\mathrm{exp}\left(\beta _N\left(\{\stackrel{}{r},\sigma \}\right)\right)`$ (7)
$`\times \delta \left(uV^1\mathrm{\Phi }(\{\stackrel{}{r},\sigma \})\right){\displaystyle \underset{\sigma }{}}\delta \left(\rho (\sigma )V^1N(\sigma )\right).`$ (8)
Integrating over the energy fluctuations yields the probability distribution function of the density distribution

$$p_V[\rho (\sigma )]=p_V[\rho (\sigma ),u]𝑑u.$$
(9)
Our specific concern is with the average form of $`\rho (\sigma )`$, given by

$$\overline{\rho }(\sigma )=\rho (\sigma )p_V[\rho (\sigma )]𝑑\rho (\sigma ).$$
(10)
Given a prescribed chemical potential distribution $`\mu (\sigma )`$ and temperature $`\beta `$, the form of $`\overline{\rho }(\sigma )`$ can be determined by simulation. Except in the ideal gas limit, however, no exact relationship between $`\overline{\rho }(\sigma )`$ and $`\mu (\sigma )`$ will generally be available. Thus one cannot (from a ‘bare’ GCE simulation), readily determine that $`\mu (\sigma )`$ corresponding to a particular desired target density distribution $`\rho _t(\sigma )`$. Subject to certain restrictions however, this can be achieved via use of histogram extrapolation.

The key idea of histogram extrapolation  is that a measured distribution $`p_V[\rho (\sigma ),u]`$ accumulated at one set of model parameters $`\mu (\sigma ),\beta `$ can be reweighted to yield estimates of the distribution appropriate to other parameters $`\mu ^{}(\sigma ),\beta ^{}`$. In its simplest form the reweighting is given by

$$p_V^{}[\rho (\sigma ),u|\mu ^{}(\sigma ),\beta ^{}]=wp_V[\rho (\sigma ),u|\mu (\sigma ),\beta ]$$
(11)
where the reweighting factor $`w`$ takes the form

$$w=\mathrm{exp}\left(\underset{i=1}{\overset{N}{}}\left[\beta ^{}\mu ^{}(\sigma _i)\beta \mu (\sigma _i)\right]V\left(\beta ^{}\beta \right)u\right)$$
(12)
By tuning the form of $`\mu ^{}(\sigma )`$ and the value of $`\beta ^{}`$ within the reweighting scheme, it is possible to scan the space of $`\overline{\rho }(\sigma )`$, thereby “homing in” on the target density distribution. To this end it is expedient to define a cost function measuring the deviation of $`\overline{\rho }(\sigma \mu ^{}(\sigma ),\beta ^{})`$ from the target form

$$\mathrm{\Delta }(\mu ^{}(\sigma ),\beta ^{})\overline{\rho }(\sigma )\rho _t(\sigma )\gamma (\sigma )𝑑\sigma .$$
(13)
Here, for numerical convenience, we have incorporated a weight function $`\gamma (\sigma )`$ into our definition, the role of which (as described below) is to ensure that comparable contributions to the cost function arise from all sampled regions of the $`\sigma `$-domain. Within this framework, the task of determining those values of $`\mu ^{}(\sigma ),\beta ^{}`$ that yield the target distribution $`\rho _t(\sigma )`$ reduces to that of functionally minimizing the cost function $`\mathrm{\Delta }`$ with respect to $`\mu ^{}(\sigma ),\beta ^{}`$. As we describe in the next section, this is achievable using standard algorithms.

### B Implementation

We have employed Monte Carlo (MC) simulation within the grand canonical ensemble to study the dilution line properties of systems of size-disperse hard spheres. This section details the principal aspects of our simulation and analysis procedure.

#### 1 Program, data structure and acquisition

Our simulated system comprises a variable number of hard spheres contained within a periodic box of volume $`V=L^d`$. The dimensionless polydisperse variable $`\sigma `$ was taken to be the particle diameter expressed in units of the mean diameter (see section III). An upper bound $`\sigma _c`$ was placed on the permitted range of diameters, and the simulation volume was partitioned into an array of $`l^3`$ cells each of linear dimension $`\sigma _c`$, so that $`L=l\sigma _c`$. This strategy aids efficient identification of particle interactions by ensuring that interacting particles occupy either the same cell or directly neighboring ones.

The grand canonical ensemble Monte Carlo (GCMC) algorithm employed has a Metropolis form  and invokes four types of operation: particle displacements, particle insertions, particle deletions and particle resizing; each is attempted with equal frequency. Specific to the polydisperse case is the resizing operation which entails attempting to change the diameter of a nominated particle by an amount drawn from a uniform random deviate constrained to lie in some prescribed range. This range (maximum diameter step-size) is chosen to provide a suitable balance between efficient sampling and a satisfactory acceptance rate at the prevailing number density. As regards the remaining types of moves, these proceed in a manner similar to the monodisperse case , except that for insertion attempts the new particle diameter is drawn with uniform probability from the range $`\sigma [0,\sigma _c]`$.

As primary input, the program takes the chemical potential distribution $`\mu (\sigma )`$, which is required for the accept/reject Monte Carlo lottery. This distribution is stored in the form of a histogram, constructed by dividing the truncated interval $`0<\sigma <\sigma _c`$ into a prescribed number $`M`$ of sub-intervals or “bins” . Thus all particles whose diameters fall within the scope of a given bin are associated with the same value of the chemical potential.

The principal observable of interest is the joint probability distribution $`p_V[\rho (\sigma ),u]`$ (c.f. Eq. 8). Operationally it is impractical to construct the full form of this distribution in a simulation because to do so would entail forming a histogram over histograms–the memory storage costs of which would be prohibitive. The procedure adopted, therefore, was to sample the instantaneous density histogram $`\rho (\sigma )`$ (discretized in the same manner as $`\mu (\sigma )`$), and append successive measurements of this quantity to a file . The set of all such samples constitutes a list representing a sequential history of the individual data measurements . This data list is post-processed by an analysis program which reads in each of the individual list entries and averages over all of them to construct a histogram approximation to the average density distribution $`\overline{\rho }(\sigma )`$. If desired, the analysis program additionally implements histogram reweighting of the data in order to enable extrapolation to neighboring values of the model parameters. This extrapolation is achieved by assigning a weight to each list entry of the form given by equation 12. The complete set of weights permits construction of the reweighted histogram.

#### 2 Tracking a phase space path

The strategy we have adopted is to traverse the phase space path of interest in a stepwise fashion, utilizing histogram extrapolation to proceed from one step to the next. In the general case of particles having a finite potential, the phase space path may involve changes in the temperature as well the form of $`\rho (\sigma )`$. For simplicity of illustration, however, let us presume that the path is isothermal (i.e $`\beta =\mathrm{constant}`$) and that the form of the chemical potential distribution, $`\mu ^{(0)}(\sigma )`$ say, is known at some arbitrary point $`\rho ^{(0)}(\sigma )`$ along the path. The procedure is then as follows. From simulations at the known state point, data for $`p_V[\rho (\sigma ),u]`$ can be accumulated directly. Histogram reweighting is then applied to this data to extrapolate some distance along the path to a new point $`\rho ^{(1)}(\sigma )`$, and to provide an estimate of the corresponding form of $`\mu ^{(1)}(\sigma )`$. The latter quantity is then employed in a fresh simulation, the results of which are extrapolated to a point further along the path, and so on. By iterating this procedure $`\rho ^{(i)}(\sigma )\rho ^{(i+1)}(\sigma )`$, $`\mu ^{(i)}(\sigma )\mu ^{(i+1)}(\sigma )`$, one traces out the entire phase space path.

The implementation of the extrapolation stage necessitates a prior choice for the step size, that is the difference between the measured $`\rho (\sigma )`$ and the next target $`\rho _t(\sigma )`$. The magnitude of this difference should be chosen to be as large as possible, consistent with remaining within the range of reliable extrapolation. A good indicator that this is in fact the case is that the individual densities of the target distribution $`\rho _t(\sigma )`$ each overlap with the range of typical fluctuations appearing in the simulation data $`\rho (\sigma )`$ .

Once a suitable step size has been determined, the extrapolation procedure proceeds by minimizing (within the reweighting scheme) the cost function $`\mathrm{\Delta }`$ introduced in eq. 13. For all but the lowest densities, this task is complicated by the existence of strong coupling between the $`\mu `$ variables, deriving from the fact that the number density for each $`\sigma `$ depends on the whole chemical potential distribution. Fortunately, efficient algorithms for performing multi-dimensional functional minimization are widely available  and, at least for the cases we have considered, appear to operate effectively. The sole difficulty encountered was that, on occasion, the minimization failed to fully converge for $`\sigma `$ values in the wings of $`\rho _t(\sigma )`$. To remedy this problem, a weight function $`\gamma (\sigma )`$ was incorporated in the cost function (c.f. eq. 13), the purpose of which is to enhance the contribution to $`\mathrm{\Delta }`$ from $`\sigma `$ values for which $`\rho _t(\sigma )`$ is small. Good results were obtained by setting $`\gamma (\sigma )[\rho _t(\sigma )]^1`$ .

It should be stressed that the methodology set out above presumes the availability of a form of $`\mu (\sigma )`$ for some starting point on the phase space path of interest. This may be obtained straightforwardly if the path passes through a region of low density where reliable analytical estimates for $`\mu (\sigma )`$ can be employed. Otherwise the method must be bootstrapped by other means. One simple but effective approach for achieving this operates as follows. Starting from some initial guess for the desired $`\mu (\sigma )`$ (e.g. $`\mu ^0(\sigma )=\mathrm{ln}[\rho _t(\sigma )]`$) a series of short simulations are carried out in which $`\mu (\sigma )`$ is iterated according to

$$\mu ^{(m+1)}(\sigma )=\mu ^{(m)}(\sigma )+\left(\frac{\mathrm{ln}\left(\rho _t(\sigma )\right)}{\mathrm{ln}\left(\rho ^{(m)}(\sigma )\right)}\right)^\delta ,$$
(14)
where $`0<\delta <1`$ is a damping factor, the value of which may be tuned to optimize convergence. Although one could certainly envisage more sophisticated and efficient schemes, we have found this method to operate satisfactorily.

## III Dilution line studies of polydisperse hard sphere fluids

### A System and simulation details

We have obtained the dilution line properties (cf. sec. I) of size-disperse hard sphere  fluids with diameter distribution $`N(\sigma )`$ assigned one of two forms:

* Schulz distribution.
* Log-normal distribution.

These distributions are conveniently expressed in terms of a normalized size function $`n(\sigma )=N(\sigma )/N`$. For the Schulz, this takes the form

$$n(\sigma )=\frac{1}{z!}\left(\frac{z+1}{\overline{\sigma }}\right)^{z+1}\sigma ^z\mathrm{exp}\left[\left(\frac{z+1}{\overline{\sigma }}\right)\sigma \right],$$
(15)
where $`\overline{\sigma }`$ is the average particle diameter and $`z`$ is a parameter controlling the width of the distribution. For the log-normal distribution, the size function is given by

$$n(\sigma )=\frac{1+W^2}{\overline{\sigma }\sqrt{2\pi \mathrm{ln}(1+W^2)}}\mathrm{exp}\left(\frac{[\mathrm{ln}(\sigma /\overline{\sigma })+(3/2)\mathrm{ln}(1+W^2)]^2}{2\mathrm{ln}(1+W^2)}\right),$$
(16)
with $`W`$ the standard deviation in units of $`\overline{\sigma }`$. Note that both of these distribution are normalized, that is $`_0^{\mathrm{}}n(\sigma )𝑑\sigma =1`$, and vanish as $`\sigma 0`$, implying a natural lower limit to $`\sigma `$. By contrast, there is no finite upper limit and consequently, for simulation purposes, it was necessary to impose an upper bound (cut off) $`\sigma _c`$ (see also sec. II B).

We have studied three distinct size distributions—two of the Schulz form and one of the log-normal form. For the Schulz distribution, width-parameter values of $`z=15`$ and $`z=5`$ were considered, with cutoff values $`\sigma _c=3`$ and $`\sigma _c=4`$ respectively. For the log-normal distribution, the single case $`W=2.5`$ with $`\sigma _c=12`$ was studied. In each instance we set $`\overline{\sigma }=1`$. Cell array sizes (cf. sec. II B 1) of linear dimension $`l=3,4,5`$ were used. In absolute dimensionless units ($`L=l\sigma _c`$) these correspond to $`L=9,12,15`$ respectively for the $`z=15`$ Schulz, to $`L=12,16,20`$ for the $`z=5`$ Schulz, and to $`L=36,48,60`$ for the log-normal distribution. The histogram discretization parameter (cf. sec. II B) was set to $`M=75`$ and $`M=100`$ for the $`z=15`$ and $`z=5`$ Schulz distributions respectively, and to $`M=120`$ for the log-normal distribution.

The average density distribution can be expressed in terms of the normalized size function multiplied by the overall number density $`\rho _0=N/V`$, i.e.

$$\overline{\rho }(\sigma )\rho _0n(\sigma ).$$
(17)
The procedure for tracking the dilution line utilized this overall density as a measure of the location along the line. The tracking procedure was initiated at a small value of $`\rho _0`$ by approximating the chemical potential distribution according to the ideal gas relation $`\mu (\sigma )=\mathrm{ln}\overline{\rho }(\sigma )`$ . Histogram extrapolation was applied to the resulting simulation data in order to refine the initial estimate of $`\mu (\sigma )`$. Thereafter the strategy described in sec. II B 2 was implemented to follow the dilution line to higher densities. At each step this entailed setting a target form for $`\overline{\rho }(\sigma )`$ corresponding to a value of $`\rho _0`$ larger than that used at the previous step. The cost function measuring the discrepancy between $`\overline{\rho }(\sigma )`$ and the target was then minimized to yield an estimate for the appropriate $`\mu (\sigma )`$. On efficiency grounds, this minimization was performed in two stages; an initial approximation to $`\mu (\sigma )`$ was obtained from a one-dimensional minimization in which the activity distribution $`\mathrm{exp}(\mu (\sigma ))`$ was multiplied by an overall factor. Thereafter the complete functional minimization was performed to yield a more accurate form for $`\mu (\sigma )`$.

Although (within the specific context of the tracking procedure), $`\rho _0`$ provides a convenient measure of the location on the dilution line, it conveys little information regarding the degree of packing within a polydisperse system. We shall therefore find it convenient to quote values for the overall volume fraction of the system given by

$$\eta _0^{\sigma _c}𝑑\sigma \frac{\pi }{6}\sigma ^3\overline{\rho }(\sigma ).$$
(18)
It is this quantity (rather than $`\rho _0`$) which features in the presentation of our results, to which we now turn.

### B Results

Owing to the computational complexity of the simulations, we have obtained the complete dilution line for each fluid only for the $`l=3`$ cell array size. For the larger system sizes, only a few spot measurements were made along the dilution line, each of which was bootstrapped by appeal to the measured $`\mu (\sigma )`$ obtained from the $`l=3`$ systems.

The dilution lines for the three fluids were tracked to the highest computationally accessible volume fraction, terminating only once the relaxation timescale became excessively large. The maximum volume fractions attained were $`\eta =0.445`$ and $`\eta =0.40`$ for the Schulz density distributions with $`z=15`$ and $`z=5`$ respectively, and $`\eta =0.33`$ for the log-normal case. Snapshot configurations for all three fluids are shown in fig. 1 for $`\eta `$ values slightly below the maximum attained in each case. We mention in passing that, at least for the cases of the two Schulz distributions, the maximum volume fractions reached are somewhat larger than would be readily attainable for GCMC simulations of monodisperse hard spheres. The latter, of course, become highly inefficient at large $`\eta `$ because the acceptance rate for particle insertions falls rapidly as the free volume diminishes. Whilst the same is true in the polydisperse context for insertions of large particles, space can often be found for placing a small particle. This facilitates fluctuations in the overall particle number $`N`$, while resizing operations (whose efficiency at large $`\eta `$ is typically greater than that of inserting large particles) ensure continued proper sampling of the density distribution.

The measured forms for $`\mu (\sigma )`$ as a function of $`\eta `$ differ qualitatively between the Schulz and the log-normal density distributions, and accordingly we discuss them separately. Beginning with the Schulz case, figs. 2 and 3 show the measured $`\overline{\rho }(\sigma )`$ and the corresponding $`\mu (\sigma )`$ for the two fluids at a selection of volume fractions along their respective dilution lines. The range of $`\sigma `$ values shown is that for which the simulations delivered data of reasonable statistical quality. One notes that for both the $`z=15`$ and $`z=5`$ cases, the tails marking the large–$`\sigma `$ vestiges of the distributions fall considerably short of the respective cutoffs $`\sigma _c`$. Indeed, we have verified that in the course of the simulations, no particles of diameters approaching the cutoff diameter occurred on the simulation timescale, implying that our data is unaffected by its imposition. As regards the forms of the chemical potential distributions, one sees that for small $`\eta `$ they display a maximum near the peak in $`\overline{\rho }(\sigma )`$—behavior that is of course mandated at sufficiently low $`\eta `$ by the known properties of the ideal gas limit. For larger $`\eta `$, however, the peak is lost and $`\mu (\sigma )`$ increases monotonically. The increase of $`\mu (\sigma )`$ in the regime of large $`\sigma `$ indicates that the excess chemical potential (measuring the work associated with inserting a sphere) grows faster with $`\sigma `$ than the decrease in the ideal contribution associated with the decay of $`\overline{\rho }(\sigma )`$.

Turning now to the log-normal case, fig. 4(a) demonstrates that $`\overline{\rho }(\sigma )`$ decays extremely slowly with increasing $`\sigma `$. The peak in the distribution therefore occurs at much smaller $`\sigma `$ than for the two Schulz forms although the average diameter $`\overline{\sigma }`$ is identical in all three cases. As a practical consequence of the slow decay (and notwithstanding the imposition of a very large value of $`\sigma _c`$), ranges of particle diameters extending up to the cutoff value were observed in the system. Indeed it was not feasible within the computational constraints to utilize a cutoff for which this did not occur, and hence truncation effects are always significant in this system—a point to which we shall return below. The measured forms of $`\mu (\sigma )`$ for the log-normal fluid are shown in fig. 4(b). In contrast to the Schulz distributions, they display (for all accessible $`\eta `$) a narrow maximum close to the peak in $`\overline{\rho }(\sigma )`$ at small $`\rho `$. Thereafter with increasing $`\sigma `$ there is a slow fall to a broad minimum, whereafter $`\mu (\sigma )`$ increases strongly.

Having outlined the main qualitative features of the relationship between $`\overline{\rho }(\sigma )`$ and $`\mu (\sigma )`$, it is instructive to perform a detailed comparison between our measurements and the predictions of analytical equations of state (EOS) appearing in the literature. For hard spheres, two commonly used equations are that due to Boublik, Mansoori, Carnahan, Starling and Leland (BMCSL)  based on the Carnahan-Starling equation for monodisperse hard spheres, and that due to Salacuse and Stell , based on the Percus-Yevick (PY) theory. Both are reproduced in appendix A, and express $`\mu (\sigma )`$ in terms of an expansion to third order in $`\sigma `$, with coefficients given in terms of the first three moments of $`\overline{\rho }(\sigma )`$. We have compared the predictions of the EOS for each of the three fluids studied, with the finite-size simulation data at three values of $`\eta `$, namely a low, a moderate and a high value. We describe our findings for each fluid in turn.

The results for the Schulz distribution with $`z=15`$ at the low volume fraction $`\eta =0.056`$ are shown in fig. 5(a). At first sight there is good agreement between the $`L=9`$ and $`L=12`$ simulation data and both the BMCSL and PY equations of state over the entire region of $`\sigma `$. However, closer inspection reveals appreciable discrepancies between theory and simulation, not visible on the scale of $`\mu (\sigma )`$. These are apparent if one suppresses the dominant ideal gas contribution to expose the excess chemical potential, given by (see appendix A)

$$\mu _{\mathrm{ex}}(\sigma )=\mu (\sigma )\mathrm{ln}[\rho _0n(\sigma )].$$
(19)
This quantity is plotted in the inset of fig. 5(a), from which one sees that compared to the simulation results, the EOS slightly underestimate $`\mu (\sigma )`$.

Fig. 5(b) shows the results for the $`z=15`$ Schulz distribution at the moderate volume fraction $`\eta =0.257`$. Again there is good agreement between the $`L=9`$ and $`L=12`$ simulation data suggesting that finite-size effects are insignificant. Here, however, discrepancies between the BMCSL and PY equations of state are larger than at the lower value of $`\eta `$, being evident on the scale of the absolute chemical potential. One sees that both EOS significantly underestimate $`\mu (\sigma )`$ for large $`\sigma `$, although they agree quite well with one another. A similarly high level of agreement between the data from the $`L=9`$ and $`L=12`$ system sizes is manifest at the higher volume fraction $`\eta =0.426`$, fig. 5(c). Here the PY equation of state is seen to fare somewhat better than the BMCSL equation although both underestimate $`\mu (\sigma )`$ substantially towards the upper end of the $`\sigma `$ range.

A similar picture emerges for the Schulz distribution with $`z=5`$ (fig. 6). Again both equations of state underestimate $`\mu (\sigma )`$, even at the lowest volume fraction, although the PY equation corresponds significantly more closely to the simulation results at high volume fractions than does the BMCSL equation. We could again discern no evidence for appreciable finite-size effects in the simulation results.

In seeking to compare the results for the log-normal system with the EOS predictions it is essential to bear in mind the importance of truncation effects. The moments of a truncated log-normal distribution can differ dramatically from those of the full distribution even for large values of $`\sigma _c`$. In order to facilitate a fair comparison with theory, the analytic form of $`\mu (\sigma )`$ must therefore be calculated using the moments of the same truncated distribution as employed in the simulations. The results of performing this comparison are presented in fig. 7. At low volume fraction, there is good agreement between the EOS predictions and the simulation results, while at high volume fraction the EOS underestimate the measured $`\mu (\sigma )`$ considerably, with the degree of discrepancy increasing towards the tail of the distribution. Once again we could discern no evidence of finite-size effects within the statistical uncertainties of our data.

The results presented above indicate that the BMCSL and PY equations fail to provide a quantitatively accurate description of the chemical potential distribution, particularly when the volume fraction is large. It is instructive to examine the implications of this finding for calculations of densities, which depends very sensitively (indeed exponentially so) on the chemical potential. To this end we have studied the degree to which the form of $`\mu (\sigma )`$ calculated via the BMCSL equation from some prescribed $`\rho (\sigma )`$, actually yields this density distribution when input to a simulation. The results of performing this comparison are shown in figures 8(a) – 8(c) for each of the three fluids at a high volume fraction. In each instance, the solid line shows the input density distribution $`\rho (\sigma )`$ from which $`\mu (\sigma )`$ is calculated. The data points show the simulation results obtained using this form of $`\mu (\sigma )`$ for $`3`$ different system sizes. As the figures clearly demonstrate, the measured form of $`\rho (\sigma )`$ deviate substantially from the prediction.

Finally in this section, we examine the moment structure of the excess chemical potential $`\mu _{\mathrm{ex}}(\sigma )`$ given by eq. 19. Both the BMCSL and PY equations assume that $`\mu _{\mathrm{ex}}(\sigma )`$ is expressible in terms of a cubic polynomial in $`\sigma `$. Bartlett has reached the same conclusion using geometrical arguments inspired by scaled particle theory . We have investigated this proposal by fitting our data to the expression

$$\mu _{\mathrm{ex}}(\sigma )=\mathrm{ln}(1\eta )+\alpha _1\sigma +\alpha _2\sigma ^2+\alpha _3\sigma ^3,$$
(20)
where the constant term is fixed by the requirement that in the limit $`\sigma 0`$ the probability of inserting a sphere is proportional to $`1\eta `$. We find that all our data are fitted very well by this expression; fig. 9(a) shows a typical fit for the case of the Schulz fluid with $`z=5`$ at $`\eta =0.377`$. Also shown, (fig. 9(b)) are the fit coefficients $`\alpha _1,\alpha _2,\alpha _3`$ for various values of $`\eta `$, together with the predictions of the BMCSL equation of state. One sees that at high volume fraction, the BMCSL underestimates all coefficients, the relative discrepancy being marginally larger for the $`\alpha _1`$ coefficient than for the others.

## IV Discussion and conclusions

In summary, we have presented a grand canonical simulation method for studying polydisperse fluids. The method utilizes histogram extrapolation techniques to track efficiently an arbitrary path of interest through the space of the density distribution $`\rho (\sigma )`$. We have applied it to the specific problem of obtaining the dilution line properties of size-disperse hard spheres. It should also prove useful in studying more general species of polydisperse fluids and their phase transitions, both in the bulk and confined geometries. We intend to report on such extensions in future communications.

Previous simulation studies of polydisperse fluids have generally operated within a semi-grand canonical framework  in which a fixed number of particles are studied either at constant pressure (see eg. ref. ) or within a Gibbs Ensemble MC scheme . In common with the present work, these studies utilized a fluctuating particle size distribution, realized by means of MC resizing moves controlled by a chemical potential distribution. In contrast to our approach, however, thermodynamic properties were studied as a function of the shape of the activity distribution $`\mathrm{exp}(\mu (\sigma ))`$; no constraints were placed on the conjugate density distribution, which was consequently free to adopt whichever functional form minimized the free energy for the imposed $`\mu (\sigma )`$. In view of this, one might question the extent to which the simulation results reflect the actual situation in realistic systems. One situation in which the lack of a constrained density distribution might be relevant is the interesting issue of the influence of polydispersity on the freezing of hard spheres . This was investigated in ref. , using semi-grand canonical MC and Gibbs-Duhem integration. It would be of considerable interest to harness the approach described in the present work to investigate the effects on the freezing behaviour of varying the width of the density distribution whilst constraining its shape to some physically realistic form.

In the present study, attention was focused on fluids having a high degree of polydispersity. The motivation for this choice was two-fold. First, models exhibiting a wide density distribution provide a suitably testing challenge against which to assess the effectiveness of our method. Second, there exist in the literature a number of interesting predictions concerning the role of depletion forces in highly polydisperse systems. For instance, it has been suggested by several authors that attractive depletion forces might engender novel phase transitions in polydisperse hard spheres . Specifically, Sear has suggested that a fluid of hard spheres having a log-normal size distribution will be unstable with respect to crystallization of the large particles at all finite volume fractions . By contrast, Cuesta  has predicted that a sufficiently wide log-normal distribution will exhibit fluid-fluid phase separation at some finite density. It seems likely that for a truncated size distribution a phase transition will not occur for arbitrarily small volume fraction because of the absence of the largest particles which mediate the greatest depletion forces. Thus our simulation results for the truncated distribution are able neither to conclusively confirm nor refute these predictions. It suffices to say that we observed no evidence of crystallization in the particular truncated log-normal fluid studied, up to a volume fraction of $`\eta =0.33`$. Similarly, within the range of accessible volume fractions, no evidence for phase transitions was observed in either of the two Schulz fluids studied.

Turning finally to the comparison between our simulation results and the predictions of the polydisperse equations of state, we find that neither the BMCSL nor PY equations offer a quantitatively accurate description of the thermodynamics of hard spheres for large polydispersity and at high volume fraction . Both equations underestimate $`\mu (\sigma )`$ at all fluid densities over the entire range of $`\sigma `$, implying that they underestimate the Helmholtz free energy density $`f=_0^{\rho _0}\mu (\sigma |\rho _0^{}n(\sigma ))𝑑\sigma 𝑑\rho _0^{}`$ and hence overestimate the stability of the fluid. The magnitude of this overestimate becomes more pronounced the greater the volume fraction. Interestingly we find that in this regime the PY equation performs appreciably better than the BMCSL equation despite the fact (cf. appendix A) that the latter derives from a monodisperse hard sphere equation of state which has been found to be superior to the PY approximation. It remains to be seen to what extent the overestimate of fluid stability impinges on the results of existing calculations of depletion-force induced phase separation based on the BMCSL and PY approximations . In any case, our results should provide a useful testing ground for any future improvements to the existing polydisperse equations of state .

###### Acknowledgements.

NBW acknowledges helpful discussions with Marcus Müller. This work was supported by The University of Liverpool Research Development Fund and the UK CCP5 group. PS acknowledges financial support through EPSRC grant GR/R52121.

## A Polydisperse hard sphere equations of state

The equations of state that we quote here express the chemical potential as a function the sphere diameter $`\sigma `$. The expression due to Salacuse and Stell  is a generalization of the Percus-Yevick result for monodisperse hard spheres:

$`\mu _{\mathrm{PY}}(\sigma )`$ $`=`$ $`\mathrm{ln}[\rho _0n(\sigma )]\mathrm{ln}(1\zeta ^3)+{\displaystyle \frac{\sigma ^3\zeta _0+3\sigma ^2\zeta _1+3\sigma \zeta _2}{1\zeta _3}}`$ (A1)
$`+{\displaystyle \frac{3\sigma ^3\zeta _1\zeta _2+9\sigma ^2\zeta _2^2/2}{(1\zeta _3)^2}}+{\displaystyle \frac{3\sigma ^3\zeta _2^3}{(1\zeta _3)^3}}`$ (A2)
The BMCSL equation of state , on the other hand, generalizes the Carnahan-Starling expression, which for monodisperse hard spheres is more accurate than the Percus-Yevick result. It is given by

$`\mu _{\mathrm{BMCSL}}(\sigma )`$ $`=`$ $`\mathrm{ln}[\rho _0n(\sigma )]+\left(3\sigma ^2\zeta _2^2/\zeta _3^22\sigma ^3\zeta _2^3/\zeta _3^3\right)\mathrm{ln}(1\zeta ^3)`$ (A3)
$`+{\displaystyle \frac{\sigma ^3(\zeta _0\zeta _2^3/\zeta _3^2)+3\sigma ^2\zeta _1+3\sigma \zeta _2}{1\zeta _3}}`$ (A4)
$`+{\displaystyle \frac{\sigma ^3(3\zeta _1\zeta _2\zeta _2^3/\zeta _3^2)+3\sigma ^2\zeta _2^2/\zeta _3}{(1\zeta _3)^2}}+{\displaystyle \frac{2\sigma ^3\zeta _2^3}{\zeta _3(1\zeta _3)^3}}`$ (A5)
where $`\zeta _n=\pi \rho _0m_n/6`$ and $`m_n`$ is the $`n`$-th moment of $`\overline{\rho }(\sigma )`$; note that $`\zeta _3=\eta `$, the volume fraction of hard spheres.

Both equations of state include the ideal contribution to the chemical potential, $`\mu _{\mathrm{ideal}}(\sigma )=\mathrm{ln}[\rho _0n(\sigma )]`$. For general temperature — remember that we set $`\beta =1`$ for our hard sphere system — this contribution would read $`\beta \mu _{\mathrm{ideal}}(\sigma )=\mathrm{ln}\overline{\rho }(\sigma )`$. Because there has been some discussion in the literature regarding how the ideal chemical potential for polydisperse systems should be assigned (see e.g. ), it may be helpful to note that, within the grand canonical framework, $`\mu _{\mathrm{ideal}}(\sigma )`$ is unambiguously defined via the grand canonical partition function. Indeed, for an ideal system ($`_N=0`$), one readily finds from eq. 1 that $`\overline{\rho }(\sigma )=\mathrm{exp}(\beta \mu (\sigma ))`$, giving $`\beta \mu _{\mathrm{ideal}}(\sigma )=\mathrm{ln}\overline{\rho }(\sigma )`$ as before. One may worry about dimensions in the argument of the log here; but the definition in eq. 1 already implies that dimensionless units are used for lengths and for $`\sigma `$. Were this not the case, the integrals over $`\stackrel{}{r}_i`$ and $`\sigma _i`$ would need to be normalized by a unit volume $`v_0`$ and a unit value $`\sigma _0`$ for the polydisperse attribute, and the ideal chemical potential would read $`\beta \mu _{\mathrm{ideal}}(\sigma )=\mathrm{ln}[v_0\sigma _0\overline{\rho }(\sigma )]`$, with the argument of the log now manifestly dimensionless. The normalization constant $`v_0\sigma _0`$ could in fact be made dependent on $`\sigma `$; this would just give a $`\sigma `$-dependent shift in the zero of the chemical potential scale.


format: text

text:
# Present status of IGEX dark matter search at Canfranc Underground Laboratory

## 1 Introduction

Recent cosmological observations and robust theoretical arguments require an important Dark Matter component ($`\mathrm{\Omega }_{DM}2530\%`$) in our universe, which is supposed to be made mostly of non-baryonic particles. Weakly Interacting Massive (and neutral) Particles (WIMPs), which are favourite candidates to such non-baryonic component, could fill the galactic halos accounting for the flat rotation curves which are measured for many galaxies. They could be detected by measuring the nuclear recoil produced by their elastic scattering off target nuclei in a suitable detector . In this talk new WIMP constraints derived from a germanium detector of the IGEX collaboration are presented. They improve previous limits obtained with Ge ionization detectors, and enter by the first time the so-called DAMA region (corresponding to a WIMP supposedly responsible for the annual modulation effect found in the DAMA experiment ) without using mechanisms of background rejection, but relying only on the ultra-low background achieved.

## 2 Experiment

The IGEX experiment , optimized for detecting 76Ge double-beta decay, has been described in detail elsewhere. One of the IGEX detectors of 2.2 kg (active mass $``$ 2.0 kg), enriched up to 86 % in 76Ge, is being used to look for WIMPs interacting coherently with the germanium nuclei. The Ge detector and its cryostat were fabricated following state-of-the-art ultralow background techniques and using only selected radiopure material components (see Ref ).

The detector shielding has been recently modified with respect to that of the previous set-up of Ref. . The improvements concern basically the neutron moderator. We have doubled its thickness (40 cms of polyethylene and borated water tanks) and now it much better covers the whole set-up, due to the removal of the other detectors (RG-I, RG-III and COSME), which dewars did not allow us to perfectly close the polyethylene wall. These changes were motivated by a previous study of the possible sources of background in IGEX based on several simulations, which pointed out that the neutrons from the surrounding rock could contribute considerably to the low energy background. For more details on the IGEX shielding we refer to  and .

In addition to the data acquisition system used in previous runs (described in ), a specific pulse shape analysis has been implemented for the data set presented in this talk. The pulse shapes of each event before and after amplification are recorded by two 800 MHz LeCroy 9362 digital scopes. These are analyzed one by one by means of a method based on wavelet techniques which allows us to assess the probability of this pulse to have been produced by a random fluctuation of the baseline. This probability is used as a criterium to reject events coming from electronic noise or microphonics. According to the calibration of the method, it works very efficiently for events above 4 keV.

## 3 Results and prospects

The results presented in this talk are from a recent run of an exposure of Mt=80 kg days with the modified set-up previously mentioned. The spectrum obtained is shown in Figure 1 compared with the previous IGEX published spectrum of Ref. .

The energy threshold of the detector is 4 keV and the FWHM energy resolution at the 75 keV Pb X-ray line was of 800 eV. The background rate recorded was $``$0.21, $``$0.10 and $``$0.04 c/keV/kg/day between 4–10 keV, 10–20 keV, and 25–40 keV respectively. As it can be seen, the background below 10 keV has been substantially reduced (about a factor 50%) with respect to that obtained in the previous set-up , essentially due to the improved shielding (both in lead and in polyethylene-water). This suggests that the neutrons could be an important component of the low energy background in IGEX.

The exclusion plot is derived from the recorded spectrum following the same set of hypothesis and parameters used in previous papers (see ) and is shown in Fig. 2 (thick solid line). It improves the IGEX-DM previous result (thick dashed line) as well as that of the other previous germanium ionization experiments (including the last result of Heidelberg-Moscow experiment  –thick dotted line–) for a mass range from 20 GeV to 200 GeV, which encompass that of the DAMA mass region. In particular, this new IGEX result excludes WIMP-nucleon cross-sections above 7$`\times 10^9`$ nb for masses of $``$50 GeV and enters the so-called DAMA region  where the DAMA experiment assigns a WIMP candidate to their found annual modulation signal. IGEX excludes the upper left part of this region. That is the first time that a direct search experiment without background discrimination, but with very low (raw) background, enters such region. Also shown for comparison are the contour lines of the other experiments which have entered that region.

Data collection is currently in progress and some strategies are being considered to further reduce the low energy background. If this reduction is achieved, very interesting perspectives can be set for IGEX. In Fig. 3 we plot the exclusions obtained with a flat background of 0.1 c/kg/keV/day (dot-dashed line) and of 0.04 c/kg/keV/day (solid line) down to the current 4 keV threshold for an exposure of 1 kg year. In particular, the complete DAMA region could be tested with a moderate improvement of the IGEX performances. The dashed line in Fig. 3 corresponds to a flat background of 0.002 c/kg/keV/day down to a threshold of 4 keV and 24 kg y of exposure, which are the parameters expected for GEDEON (GErmanium DEtectors in ONe cryostat), a new experimental project on WIMP detection using larger masses of natural germanium planned as an extension of the IGEX dark matter search (see ref. ). GEDEON would be massive enough  to search also for the WIMP annual modulation effect and explore positively an important part of the WIMP parameter space including the DAMA region.


text:
# Identification of the optical and quiescent counterparts to the bright X-ray transient in NGC 6440

## 1 Introduction

Currently, thirteen bright X-ray sources are known to inhabit twelve Galactic globular clusters (for a recent overview, see Sidoli et al. 2001; the 13th source was recently discovered by White & Angelini 2001). They are low-mass X-ray binaries (LMXBs) of which at least 11 harbor a neutron star. Four are transient in nature: those in Liller 1, Terzan 5 and 6, and NGC 6440 (one could argue that Terzan 1 contains a transient source as well with an activity lasting decades, see Guainazzi, Parmar & Oosterbroek 1999). Thus far, none of these transients had a securely identified optical counterpart. The prospects for finding one were best for the transient in NGC 6440 because it has a relatively low reddening of $`E(BV)=1.07`$ and a non-extreme distance of $`8.5\pm 0.4`$ kpc (Ortolani, Barbuy & Bica 1994). The transient in NGC 6440 was seen in outburst in Dec. 1971 (Markert et al. 1975; Forman, Jones & Tananbaum 1976) and Aug. 1998 (In ’t Zand et al. 1999; Verbunt et al. 2000). From ROSAT HRI localizations of two sources in quiescence and optical observations during the 1998 outburst, two candidate optical counterparts were identified, V1 and V2 (Verbunt et al. 2000). V2 was evaluated to be the most likely one with $`B22.7`$ during outburst.

Pooley et al. (2001) observed NGC 6440 for 23 ksec with the Chandra X-ray Observatory on July 4, 2000, when the transient was quiescent. They found 24 point sources with 0.5-2.5 keV luminosities between 1.8$`\times 10^{30}`$ and $`9.9\times 10^{32}`$ erg s-1, which they designate CX1 through CX24. The brightest source, CX1, coincides with the optical counterpart candidate V2.

The transient in NGC 6440 erupted again in August 2001, as detected with the All-Sky Monitor (ASM) on the Rossi X-ray Timing Explorer RXTE111 as announced by the ASM team at http://xte.mit.edu. We proceeded to secure the identification of the optical counterpart through swift, high-resolution X-ray observations, and report the results in the present paper. A preliminary report was already given by In ’t Zand et al. (2001).

## 2 Observations and source identification

The ASM (Levine et al. 1996; Bradt et al. 2000) monitors NGC 6440 for 10 months of each year at a detection threshold of about 15 mCrab in the 1.5 to 12 keV band per day of observation. The transient exceeded this threshold during two weeks starting on August 2, 2001. The peak flux was $`60\pm 10`$ mCrab. Continued ASM monitoring showed a reflaring episode starting on Sep. 10, with fluxes up to 0.2 Crab and lasting for already almost 2 months at the time of this writing. The Wide Field Cameras (WFCs) on BeppoSAX (Jager et al. 1997; In ’t Zand 2001) also observed NGC 6440 during this episode and detected at least 10 type-I X-ray bursts. As soon as the community was alerted of the outburst (on August 13), we proposed X-ray observations with Chandra and optical observations with the European Southern Observatory’s New Technology Telescope in Director’s Discretionary Time.

The Chandra (Weiskopf, O’Dell & van Speybroek 1996) target-of-opportunity observation (TOO) was approved and carried out on August 18.9 UT for an exposure time of 2.6 ksec with the Advanced CCD Imaging Spectrometer (ACIS; Garmire, Nousek & Bautz, in prep.) and no grating. Despite the anticipated heavy pile up, we chose this configuration, and not HRC, to optimize the sensitivity such that not only the transient but also a few of the other sources in NGC 6440 would be detected (ACIS has 2 to 4 times more effective area for 0.3-8 keV photons than HRC-I). Thus, the association of the transient with one of the quiescent sources would become independent from the Chandra attitude accuracy. Only a $`2\mathrm{}\times 8\mathrm{}`$ subarray of the back-illuminated S3 chip was read so that the CCD frame time was 0.8 s instead of 3.2 s and pile up effects were reduced. Figure 1 shows the obtained image, zoomed in on the globular cluster and compared to the observation during quiescence (Pooley et al. 2001). It is evident that the transient was active. Furthermore, it is closest to CX1, the brightest of the sources detected during quiescence. Based on the Chandra aspect solution for the 2001 observation (accurate to about 0.5″) and the astrometry carried out by Pooley et al. (2001) for the 2000 observation, the angular distance between the transient and CX1 is 0.14″ (the position of the transient is determined through weighing pixel position with pixel values). Of the remaining sources, CX2 and CX5 are clearly detected during the TOO as well, with photon rates in the full bandpass of $`(8\pm 2)\times 10^3`$ and $`(6\pm 1)\times 10^3`$ s-1. These rates compare well with those during the earlier observation ($`(7.6\pm 0.6)\times 10^3`$ and $`(4.3\pm 0.4)\times 10^3`$ s-1 respectively). The angular distances of the transient to CX2 and CX5 are consistent with that to CX1 in the July 2000 image within 0.28$`\pm `$0.13″. The distance is inconsistent with all of the other 23 sources from the July 2000 observation by at least 15$`\sigma `$. This confirms the transient-CX1 association.

A total of $`2.6\times 10^3`$ events was detected within 8.0″ from the source centroid. In order to estimate the flux and avoid piled-up pixels, we extracted a spectrum from the wings of the point spread function (PSF) for pixels that are between 2.5 and 8.0″ from the centroid. We ignore negligible contributions from other sources in this annulus. An absorbed power law fits the 0.7-8.0 keV spectral data well, with a photon index of $`1.2\pm 0.1`$, $`N_\mathrm{H}=(8.2\pm 1.6)\times 10^{21}`$ cm-2 and $`\chi _\nu ^2=1.40`$ ($`\nu =75`$ degrees of freedom). After a correction for the fraction of the PSF covered (9%), the flux is $`1.1\times 10^{10}`$ erg s-1cm-2 (0.5-10 keV) and the luminosity $`9\times 10^{35}`$ erg s-1 ($`1.1\times 10^{36}`$ erg s-1 unabsorbed). The statistical uncertainty is 7%, that due to PSF uncertainties an estimated 10%. The flux is roughly 10 times fainter than at the peak and $`10^3`$ times brighter than during quiescence (see below). The source shows no strong variability, with an upper limit of 20% at 90% confidence for 100 s time scales. No narrow features are apparent in the power density spectrum.

The position of CX1 (Pooley et al. 2001) is consistent with that of the candidate optical counterpart V2 in Verbunt et al. (2000). They are 0.46″ apart while the uncertainty in the optical position alone is 0.5″. We conclude that V2 is the secure optical counterpart to the transient in NGC 6440, that CX1 is the quiescent X-ray counterpart, and that the 1998 and 2001 outbursts were from the same object and not from any of the other 3 or 4 transient low-mass X-ray binaries in quiescence that are likely to be present in the cluster (Pooley et al. 2001). Most importantly, since type-I X-ray bursts were detected from the transient during both 1998 and 2001 outbursts, it is now firmly established that CX1 is a quiescent neutron star.

NGC 6440 was observed for us in Service Time with the New Technology Telescope (NTT) at La Silla 3.5 d after the Chandra observation using the Superb Seeing Imager SUSI2. The seeing was around 1″. Images were obtained in the Bessell R and B filters. The analysis of these data is complicated by proper motions of a number of stars in the field. However, the raw images do show V2 to be active again in B at a level comparable to the 1998 outburst.

## 3 The quiescent X-ray spectrum

Quiescent neutron stars (qNSs) exhibit X-ray emission at typical luminosities of $`10^{32}10^{33}`$ erg s-1 (e.g., Asai et al. 1998, McClintock 2001). A plausible explanation for this emission is thermal radiation from the atmosphere of a neutron star that is cooling from the previous outburst (Brown, Bildsten & Rutledge 1998). We extracted the spectrum of the transient during quiescence by accumulating all 251 photons within 1.5″ from the source centroid and applying background corrections as found elsewhere in the image (amounting merely to a 2.3 photons subtraction). To ensure applicability of the $`\chi ^2`$ statistic, we rebinned the data so that each channel contains at least 20 photons except for the highest-energy bin which contains 15 photons. We used CIAO 2.1.3 Science Threads for the background subtraction and response matrix generation. The spectrum was fitted between 0.3 and 6.7 keV (only 2 photons were detected above this range) with a model combining a power law and a neutron star hydrogen atmosphere model as calculated for the non-magnetic case by Zavlin, Pavlov & Shibanov (1996). Standard absorption by cold interstellar matter of cosmological abundances was applied following Morrison & McCammon (1983). Due to the limited statistical quality, the fit errors for the atmosphere temperature and radius are reasonably constrained only if the mass of the neutron star, the distance, the interstellar hydrogen column density $`N_\mathrm{H}`$ and the power law photon index are fixed. We did so at values of 1.4 M (for a ’canonical’ neutron star), 8.5 kpc (this is well determined), $`N_\mathrm{H}=8.2\times 10^{21}`$ cm-2 and 1.44 (these two values were determined from the 1998 outburst; In ’t Zand et al. 1999), respectively. $`N_\mathrm{H}`$ has an uncertainty of roughly 40%. The only three parameters left free are the neutron star radius, the atmosphere temperature and the power law normalization. The resulting fit is good, with $`\chi _\nu ^2=0.758`$ (for $`\nu =9`$ degrees of freedom). The neutron star parameters are strongly dependent. The fitted temperature is $`92_{33}^{+172}`$ eV and the radius $`17_{12}^{+31}`$ km (90% confidence). The radius to an observer at infinity is $`19_8^{+29}`$ km. A graph of the spectrum is presented in figure 2 while in figure 3 the 90% confidence region in temperature/radius space is presented. If $`N_\mathrm{H}`$ is left free, it converges to $`8.2_{1.3}^{+2.0}\times 10^{21}`$ cm-2 and the range of allowed NS radii surpasses the 100 km mark. This upper radius margin is also quite sensitive to the power law index which is poorly constrained by the data.

The inferred radius is consistent with that expected for a neutron star (e.g., Lattimer & Prakash 2001). We note that the model implies an intrinsic radius that is 20% larger if the NS mass is 1.8 instead of 1.4 M. Fits with some other models are nearly as good (table 1), so the data are not discriminate in that sense. The absorbed 0.5-7 keV flux is $`(5.9\pm 0.4)\times 10^{14}`$ erg s-1cm-2, the unabsorbed flux $`(2.0\pm 0.6)\times 10^{13}`$ erg s-1cm-2 (the large error reflects the uncertainty in $`N_\mathrm{H}`$). The power law contribution in the latter flux is 20$`\pm `$6% (90% confidence); this is only 8% in 0.5 to 2 keV. The bolometric unabsorbed luminosity of the NS atmosphere is $`(2.1\pm 0.8)\times 10^{33}`$ erg s-1 (the contribution from the uncertainty in the distance is small).

## 4 Discussion

Pooley et al. (2001) find that there are 4 or 5 likely quiescent LMXB transients in NGC 6440. Similar conclusions are currently being reached from Chandra observations of other globular clusters such as 47 Tuc (Grindlay et al. 2001a), Liller 1 (Homer et al. 2001), NGC 6652 (Heinke, Edmonds & Grindlay 2001), $`\omega `$ Cen (Rutledge et al. 2001) and NGC 6397 (Grindlay et al. 2001b). Furthermore, for the first time two bright low-mass X-ray binaries have been identified in a single cluster (M15; White & Angelini 2001). Our observations show unequivocally that the second and third transient outburst in NGC 6440 (in 1998 and 2001) are from the same object. Due to a limited positional accuracy of the source of the 1971 outburst, it is not clear whether the same source was responsible for the 1971 and the 1998/2001 outbursts. At first, indirect evidence suggested that the source is different, because the duration of 1.5 to 3 months and the peak flux of 250 mCrab deviate substantially from the values for the 1998 outburst and for the initial phase of the 2001 outburst. However, at the time of this writing the 2001 outburst shows activity for already 3 months and a peak flux in excess of about 0.2 Crab units, according to ASM and WFC measurements.

The recurrence time between the two most recent outbursts is unexpectedly short. This may mean that the source has become more active, or alternatively that earlier outbursts have been missed. It is not uncommon for a previously quiet source, to become active in a series of transient outbursts (e.g. GRO J1744-28 and GRO J1655-40). However, the 1998 outburst of the transient in NGC6440 was so faint and short, that similar outbursts may easily have passed undetected. We note that the RXTE ASM and BeppoSAX WFCs became operational only in 1996. Similar coverage as the ASM in the future should resolve this issue.

This is the first time that an optical counterpart has firmly been established for an X-ray transient in a globular cluster. Furthermore, it is the first time that a quiescent neutron star has been identified without source confusion in a globular cluster. These identifications are important because they present the case of a quiescent neutron star transient in the optical and X-rays for which an accurate distance is known.

The spectral study of quiescent neutron stars is foremost interesting, as Brown et al. (1998) pointed out, because it may result in an accurate enough neutron star radius determination to be able to constrain the equation of state of matter in the star’s core. Measuring the radius to an accuracy of about 0.5 km can exclude 50% of the proposed equations of states (Lattimer & Prakash 2001). For neutron stars in globular clusters as compared to those in the Galactic disk, one source of uncertainty is almost eliminated — the distance. The distance to NGC 6440 is accurate to 5%. If other uncertainties (in the spectral parameters) had been at least an order of magnitude lower, the prospects for constraining the equation of state for the transient in NGC 6440 would have been good. However, this is not the case. This qNS would have to be observed for a few million seconds with Chandra before it becomes worthwhile to consider this analysis, and that is without considering uncertainties in the spectral models.

Given a peak bolometric luminosity between $`4\times 10^{37}`$ and $`1.6\times 10^{38}`$ erg s-1 (or 50 to 200 mCrab in 2-10 keV given the broad-band spectrum measured by In ’t Zand et al. 1999), an e-folding decay time between 5 d (Verbunt et al. 2000) and 30 d, and assuming that during accretion 100% of the gravitational potential energy is released through radiation, the recurrence time of 2.95 yr implies a time-averaged accretion rate of $`<\dot{M}>2\times 10^{11}`$ to $`4\times 10^{10}M_{}`$ yr-1. Following Brown et al. (1998), the quiescent time-averaged luminosity due to thermal emission from the cooling NS is predicted to be between $`1\times 10^{33}`$ and $`2\times 10^{34}`$ erg s-1. Our measurement of $`2.1\times 10^{33}`$ erg s-1 is consistent with this rough prediction, also if less than 100% of the liberated potential energy is released as radiation.

We thank Harvey Tananbaum for the willingness to schedule the Chandra observation at short notice, Yousaf Butt for help in planning the Chandra observation, Olivier Hainaut for carrying out the optical observations, Rob van der Meer for help in the ACIS data analysis, Slava Zavlin and George Pavlov for providing their neutron star atmosphere models in electronic form, and Jaap Schuurmans for software support. JZ acknowledges financial support from the Netherlands Organization for Scientific Research (NWO), MHvK through a fellowship of the Royal Netherlands Academy of Sciences (KNAW), DP through a National Science Foundation Graduate Fellowship, WL from NASA, and RW from NASA through a Chandra fellowship.


text:
# FORMATION OF PLANETARY NEBULA LOBES BY JETS

## 1 INTRODUCTION

Planetary nebulae (PNe) and proto-PNe possess a rich spectrum of different structures, with spherical and bipolar PNe at the two extremes of that spectrum. Bipolar PNe are defined as extreme aspherical PNe having two lobes with an ‘equatorial’ waist between them (Schwarz, Corradi & Stanghellini 1992). Several models have been suggested to explain the formation of lobes. One of the popular models a decade ago was the interacting wind model, where a fast wind blown by the central star of the PN is collimated by a previously ejected AGB stellar dense equatorial wind (e.g., Balick 1987; Soker & Livio 1989; Frank & Mellema 1994; Mellema & Frank 1995; Mellema 1995). When the equatorial to polar density ratio is very high, a bipolar nebula was supposed to be formed. However, the numerical simulations did not form bipolar PNe with very narrow waists, but rather formed elliptical PNe, or bipolar PNe with wide waists. Another problem of the interacting wind model is the high momentum and kinetic energy observed in several bipolar PNe and proto-PNe (Bujarrabal et al. 2001; see Balick 2000 for some other problems). Therefore, the simple interacting wind model cannot explain many of the observed bipolar PNe. It seems that in order to form many bipolar PNe, e.g., those with very narrow waists, a collimated fast wind (CFW) is required (Morris 1987; Soker & Rappaport 2000, hereafter SR00). Such preliminary simulations, but where the CFW was blown by the primary star, were performed by Frank, Ryu, & Davidson (1998), who indeed got very narrow waists in their simulated nebulae. Other types of models are based on magnetic fields playing a dynamical role in the mass loss process from AGB stars (e.g., García-Segura 1997; García-Segura & López 2000; Blackman et al. 2001). There are some severe problems with dynamical magnetic models (Soker & Zoabi 2002). The main problem is that the AGB has to be spun-up by a stellar companion in order to possess the required activity; hence we are in the regime of binary models.

As reviewed recently by SR00, the most promising mechanism for the formation of bipolar PNe with narrow waists is based on an accreting compact companion to the mass losing AGB star. In most cases, the companion is a white dwarf or a main sequence star (SR00). When a compact binary companion accretes matter with sufficiently high specific angular momentum, an accretion disk is formed. When the accretion rate is above some threshold (SR00), two jets are expected to be blown by the companion. If the jets are not well collimated, the outflow from the companion is termed CFW. In addition, the companion may lead to a concentrated equatorial flow (Mastrodemos & Morris 1999), further increasing the equatorial to polar density ratio. The formation of bubbles (or lobes) by jets have been studied in a variety of astrophysical flows, e.g, young stellar objects (e.g., Masson & Chernin 1993), clusters of galaxies (e.g., Reynolds, Heinz & Begelman 2001), and gamma ray burst environments (Ayal & Piran 2001). However, the condition of jets blown into AGB winds are different from those in these cases.

To study the interaction of the two winds blown simultaneously, the slow wind by the AGB star and the CFW (or two jets) by the accreting compact companion (see figs. 1 and 2 in SR00), a 3D gasdynamical numerical code is required. Such numerical simulations require heavy computer resources, hence can cover only a small number of cases. Therefore, it will be helpful to derive simple expressions than can both describe the general structure of lobes obtained by different jets, and guide future numerical simulations. This paper is aimed at reaching these objectives, by exploring part of the parameter space of jets interacting with spherical AGB winds, and under simplifying assumptions. To manifest analytical solutions, I consider two extreme cases, which are described in $`\mathrm{\S }2`$, and are studied in $`\mathrm{\S }3`$ and $`\mathrm{\S }4`$. A short summary is in $`\mathrm{\S }5`$.

## 2 JET PROPAGATION

In the present analytical study, I treat the flow at large distances from the binary system, neglecting the deflection of the jets (or CFW; see SR00). I assume an axisymmetric flow in which a jet expands along the symmetry axis and into a previously ejected spherically symmetric slow AGB wind. The slow wind density as a function of distance from the binary system is $`\rho _s=\dot{M}_s/(4\pi v_sr^2)`$, where $`v_s`$ is the slow wind speed, and $`\dot{M}_s`$ is the mass loss rate, both are constant, and where mass loss rate is defined positively. The gas inside the undisturbed jets, one at each side of the equatorial plane, expands at a constant speed $`v_j`$, the mass loss rate into each jet is $`\dot{M}_j`$, and each jet expands into a solid angle of $`\mathrm{\Omega }_j=4\pi \beta `$. The density inside an undisturbed radially expanding jet, i.e., its cross section increases as $`r^2`$, is $`\rho _j=\dot{M}_j/(4\pi \beta v_jr^2)`$. The velocity of the jet head, $`v_h`$, is determined by equating the slow wind pressure on the jet head with that of the jet material. Since thermal pressure can be neglected in the pre-shocked media, the expression is $`\rho _j(v_jv_h)^2=\rho _s(v_hv_s)^2`$. In the present analytical study I treat two extreme cases: slow-propagating and fast-propagating jets. In the slow-propagating jet case the jet head proceeds at a speed $`v_hv_j`$, while in the fast-propagating jet case the jet head does not slow much, $`v_sv_hv_j`$.

The velocity of a slow-propagating jet to be treated in the next section is given by

$`v_hv_s\left[1+\left({\displaystyle \frac{\dot{M}_jv_j}{\beta \dot{M}_sv_s}}\right)^{1/2}\right];\mathrm{for}\rho _s\rho _j.`$ (1)
I therefore assume that $`\dot{M}_jv_j\mathrm{few}\times 10\beta \dot{M}_sv_s`$, so that $`v_h210v_s`$. Since I neglect the bending of the jet by the AGB slow wind close to the binary system, the momentum flux of the jet must be (SR00) $`\dot{M}_jv_j\dot{M}_sv_s\mathrm{tan}\alpha `$, where $`\alpha `$ is the half opening angle of the jet (SR00). For $`\mathrm{tan}\alpha 1`$, $`\mathrm{tan}\alpha 2(\beta )^{1/2}`$. For $`\alpha =10^{}`$, for example, $`\beta =0.0076`$, $`2\beta ^{1/2}=0.174`$, and $`\mathrm{tan}\alpha =0.176`$. Over all, the treatment of the slow-propagating jet is applicable to the case $`(\dot{M}_jv_j)/(\dot{M}_sv_s)0.10.5`$, for a half opening angle of $`\alpha 530^{}`$. However, since, as we’ll see later, the derived expressions are not so sensitive to the different variables, the results will be applicable to a much larger parameter space.

From the expressions developed above, we can find the density of the jet shocked material, and from that, the radiative cooling time, $`t_{\mathrm{cool}}=(5/2)nkT/(n_en_p\mathrm{\Lambda })`$, where $`n`$, $`n_e`$, and $`n_p`$, are the total number, electron, and proton, post-shock densities, respectively, taken to be four times the pre-shock densities, and for the considered temperature range, which is determined by the speed of the jet material, $`\mathrm{\Lambda }10^{22}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1`$. The cooling time of the shocked jet material is obtained by substituting typical values

$`t_{\mathrm{cool}}(\mathrm{jet})17\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^3\left({\displaystyle \frac{\dot{M}_j}{10^7M_{}\mathrm{yr}^1}}\right)^1\left({\displaystyle \frac{\beta }{0.01}}\right)\left({\displaystyle \frac{r}{10^{16}\mathrm{cm}}}\right)^2\mathrm{yr}.`$ (2)
This time scale is much longer than the flow time along the jet diameter. However, the cooling time of the post shocked material should be treated differently and compared with another flow time scale. This is done for the two flow cases in the next two sections.

## 3 SLOW-PROPAGATING JETS

### 3.1 Jet Interaction with the Slow Wind

The shocked jet material will expand and form a hot low density bubble that will accelerate the dense slow wind’s material around it. The results of Castor, McCray & Weaver (1975) for a stellar wind bubble in a dense interstellar medium can be used. Using their equation (6) for the radius of the expanding bubble as a function of the expansion time of the bubble, $`t_b`$, gives

$`R_b=0.76\left({\displaystyle \frac{0.5\dot{M}_jv_j^2}{\rho _s}}\right)^{1/5}t^{3/5}=1.0\times 10^{15}\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{2/5}\left({\displaystyle \frac{\dot{M}_j}{10^7M_{}\mathrm{yr}^{1/5}}}\right)^{1/5}`$
$`\times \left({\displaystyle \frac{v_s}{15\mathrm{km}\mathrm{s}^1}}\right)^{1/5}\left({\displaystyle \frac{\dot{M}_s}{10^5M_{}\mathrm{yr}^1}}\right)^{1/5}\left({\displaystyle \frac{r}{10^{16}\mathrm{cm}}}\right)^{2/5}\left({\displaystyle \frac{t_b}{100\mathrm{yr}}}\right)^{3/5}\mathrm{cm},`$ (3)
where the local slow wind density at $`r`$ was taken in the second equality. The expansion velocity of the bubble’s surface is

$`v_b=0.6R_b/t=20\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{2/5}\left({\displaystyle \frac{\dot{M}_j}{10^7M_{}\mathrm{yr}^{1/5}}}\right)^{1/5}\left({\displaystyle \frac{v_s}{15\mathrm{km}\mathrm{s}^1}}\right)^{1/5}`$
$`\times \left({\displaystyle \frac{\dot{M}_s}{10^5M_{}\mathrm{yr}^1}}\right)^{1/5}\left({\displaystyle \frac{r}{10^{16}\mathrm{cm}}}\right)^{2/5}\left({\displaystyle \frac{t_b}{100\mathrm{yr}}}\right)^{2/5}\mathrm{km}\mathrm{s}^1.`$ (4)
The expansion time of the hot bubble, $`t_b`$, should be compared with the cooling time of the hot gas inside the bubble. The density inside the bubble is given by $`\rho _b=\dot{M}_jt_b/(4\pi R_b^3/3)`$, from which the cooling time is derived

$`t_{\mathrm{cool}}(\mathrm{bubble})30\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^2\left({\displaystyle \frac{\dot{M}_j}{10^7M_{}\mathrm{yr}^1}}\right)^1\left({\displaystyle \frac{R_b}{5\times 10^{15}\mathrm{cm}}}\right)^3\left({\displaystyle \frac{t_b}{100\mathrm{yr}}}\right)^1\mathrm{yr}.`$ (5)
The condition for the hot gas (the shocked jet material) to inflate a bubble is that the cooling time is longer than the expansion time, i.e., $`t_{\mathrm{cool}}(\mathrm{bubble})>t_b`$. Substituting for $`t_{\mathrm{cool}}(\mathrm{bubble})`$ from equation (5), and eliminating $`R_b`$ by equation (3), gives the radius outward to which the shocked jet material will not have time to cool before inflating a large bubble (which will evolve to a nebular lobe)

$`r_s4.4\times 10^{15}\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{8/3}\left({\displaystyle \frac{\dot{M}_j}{10^7M_{}\mathrm{yr}^1}}\right)^{1/3}\left({\displaystyle \frac{v_s}{15\mathrm{km}\mathrm{s}^1}}\right)^{1/2}`$
$`\times \left({\displaystyle \frac{\dot{M}_s}{10^5M_{}\mathrm{yr}^1}}\right)^{1/2}\left({\displaystyle \frac{t_b}{100\mathrm{yr}}}\right)^{1/6}\mathrm{cm}.`$ (6)
The scaling in the last equation shows that the adiabatic phase, where the cooling time of the shocked material is very long, starts at $`r10^{16}\mathrm{cm}`$. However, the parameters change from one system to another. For example, in the proto-bipolar PN OH 231.8+4.2, the mass loss rate is very high, $`2\times 10^4M_{}\mathrm{yr}^1`$ (Alcolea et al. 2002 and more references therein). If we take the mass loss rate in both the AGB wind and the jet to be 10 times higher than the values used in the scaling of equation (6), we find the adiabatic phase to start at $`r_s5\times 10^{16}\mathrm{cm}`$ from the central source. If the companion accretes via a Roche lobe overflow, then the mass loss rate into the jets can be much higher. In OH 231.8+4.2, the mass loss rate into each of the two jets could have been $`10^4M_{}\mathrm{yr}^1`$ for a short time during a Roche lobe overflow (Soker 2002). If we assume that the opening angles of the two jets were large, so the jets’ density is lower than the slow wind density, we can use equation (6). Substituting $`\dot{M}_s=10^4M_{}\mathrm{yr}^1`$, $`\dot{M}_j=10^4M_{}\mathrm{yr}^1`$, and $`v_j=500\mathrm{km}\mathrm{s}^1`$, we find from equation (6) $`r_s5\times 10^{16}\mathrm{cm}`$. This is a substantial fraction of the lobe size (Kastner & Weintraub 1995). Therefore, during a large fraction of the evolution (since the jet phase was short; Soker 2002) the flow was momentum conserving (radiative), i.e., the shocked material was cooling on a relatively short time. Indeed, Bujarrabal et al. (2002), ruled out the adiabatic case for OH 231.8+4.2, mainly on the ground of the shape of the lobes, which are not inflated as is expected in the adiabatic flow.

### 3.2 The Bubble Shape

From equation (3) it turns out that the radius of a spherical bubble increases with $`r`$ much slower than $`r`$ does. For example, at a distance from the center 32 times that was is used in equation (3) the bubble radius will be only 4 times larger, keeping the time and the other parameters unchanged. However, the bubble will not be spherical since the jet propagates outward, and a more accurate treatment is required. I assume that the bubble expands perpendicular to the symmetry axis, along the $`z`$ direction, i.e., parallel to the equatorial plane. I also take the density in planes parallel to the equatorial plane to be the same as the density where the plane intersects the symmetry axis. The distance of the plane from the equatorial plane is the coordinate $`x`$. The relevant quantity is the energy deposited per unit length along the jets propagation, which is $`e_0=0.5\dot{M}_jv_j^2/v_h`$. I use the relation between pressure and energy and the momentum equations as in Castor et al. (1975), and neglect the thermal energy of the swept slow wind material. The last assumption is adequate since I neglect the slow wind radial expansion and the dependence of density on the distance from the symmetry axis, while the assumptions allow a simple expression for the energy equation. The energy equation reads $`e=e_0e_k`$, where $`e`$ is the thermal energy per unit length (along the symmetry axis) in the bubble, and $`e_k`$ is the kinetic energy per unit length of the swept up slow wind gas. Neglecting the density variation perpendicular to the symmetry axis and the radial expansion of the slow wind gives the distance of the bubble surface from the symmetry axis

$`z_b(x)=\left({\displaystyle \frac{16}{5\pi }}{\displaystyle \frac{e_0}{\rho _s(x)}}\right)^{1/4}t^{1/2}\left({\displaystyle \frac{e_0}{\rho _s(x)}}\right)^{1/4}t^{1/2},`$ (7)
where $`x`$ and $`z`$ are the coordinates parallel and perpendicular to the symmetry axis, respectively. Substituting in the typical values used here, with $`v_h`$ by equation (1), gives

$`{\displaystyle \frac{z_b}{x}}1.3\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{1/2}\left({\displaystyle \frac{\dot{M}_j}{10^7M_{}\mathrm{yr}^1}}\right)^{1/4}\left({\displaystyle \frac{v_h}{4v_s}}\right)^{1/4}\left({\displaystyle \frac{\dot{M}_s}{10^5M_{}\mathrm{yr}^1}}\right)^{1/4}`$
$`\times \left[{\displaystyle \frac{t_b(x)}{100\mathrm{yr}}}\right]^{1/2}\left({\displaystyle \frac{x}{10^{16}\mathrm{cm}}}\right)^{1/2}.`$ (8)
For the same parameters, the velocity is

$`v_b{\displaystyle \frac{dz_b}{dt}}0.5{\displaystyle \frac{z_b}{t}}20\left[{\displaystyle \frac{t_b(x)}{100\mathrm{yr}}}\right]^{1/2}\left({\displaystyle \frac{x}{10^{16}\mathrm{cm}}}\right)^{1/2}\mathrm{km}\mathrm{s}^1.`$ (9)
Because of the assumption that density does not depend on the $`z`$ coordinate, the last two equations are less accurate for $`z_b>r`$, i.e., close to the central binary system.

It should be noted that in the last two equations the time $`t_b(x)`$ is counted from the moment the jet reaches the point $`x`$, hence $`t_b`$ is longer for lower values of $`x`$, i.e., closer to the center of the nebula. If the nebula is observed at a time $`t_{\mathrm{obs}}`$ after the jets are launched, and the jets proceed at velocity $`v_h`$, then $`t_b(x)=t_{\mathrm{obs}}x/v_h`$. Substituting this in equation (8), keeping the other variables the same, gives

$`z_b1.3\times 10^{16}\left[{\displaystyle \frac{t_{\mathrm{obs}}(x/v_h)}{100\mathrm{yr}}}\right]^{1/2}\left({\displaystyle \frac{x}{10^{16}\mathrm{cm}}}\right)^{1/2}.`$ (10)
The maximum width of the lobe, $`Z_{b\mathrm{max}}`$, occurs at

$`x(z_{b\mathrm{max}}){\displaystyle \frac{v_ht_{\mathrm{obs}}}{2}}=1.6\times 10^{17}\mathrm{cm}\left({\displaystyle \frac{\mathrm{t}_{\mathrm{obs}}}{1000\mathrm{yr}}}\right)\left({\displaystyle \frac{v_h}{100\mathrm{km}\mathrm{s}^1}}\right)\mathrm{cm},`$ (11)
assuming the jet life time was $`>t_{\mathrm{obs}}/2`$. For this time of observation, $`v_h100\mathrm{km}\mathrm{s}^1`$, and with the other values in equation (8), the maximum width is $`z/x0.7`$. Due to the decrease of the density with distance from the symmetry axis and the radial expansion of the slow wind, the lobe will actually be somewhat wider. These values for the location of the widest lobal point and its width there are typical for many young bipolar PNe. The derivation above shows that the bubble, or lobe, is narrow close to the central system, then becomes wider, and then narrows again. If the density of the material on the outer portion is too low to be observed, the observed lobe surface will only widen with distance from the equatorial plane.

There are several strong assumptions in the derivation above; in particular, $`(i)`$ neglecting the variation of the slow wind density perpendicular to the symmetry axis; $`(ii)`$ treating the slow wind as a static medium rather than a radially expanding medium; and $`(iii)`$ the assumption that the bubble surface only expands perpendicular to the symmetry axis (the $`z`$ direction). The last assumption is reasonable at the intermediate values of $`x`$, but not close to the equatorial plane (small values of $`x`$). The expansion there is more appropriately described by equation (3) and (4). Also, from equation (6) it is clear that the center of each bubble will be at a large distance from the equator. Hence, the expanding bubble surface facing the equatorial plane will push material toward the equatorial plane, substantially increasing the density there. This supports the claim of SR00 that jets (or a CFW) can compress the gas near the equatorial plane, contrary to some claims that the dense equatorial gas collimates the jets. If the bubble surface reaches the equatorial plane, it will collide with the bubble on the other side of the equatorial plane, and the lobe sides will be observed to emerge from the equatorial plane at a large distance from the center. Such a structure is observed, for example, in the PN He 2-104 (Sahai & Trauger 1998). The inflation of the bubble away from the central system implies that the bubble boundary may not reach the equatorial plane. This effect is stronger if we consider the outward motion of the slow wind, which was ignored in the derivation above. However, the distance of the bubble surface from the equatorial plane will be very small, since the bubble is starting to be inflated by the jet (or CFW) at $`r400\mathrm{AU}`$ (eq. 6). To observe a bubble not touching the equatorial plane, the nebula must be observed almost edge on. These two constraints imply that only a very small number of such bipolar PNe should be observed. Taking into account its inclination, it seem that the lobes of the bipolar-PN Hb 12 do not cross the equatorial plane (Sahai & Trauger 1998). At larger values of $`x`$ the jets may cease to exist, and a different treatment is required. Despite these approximations, the treatment above gives the basic shape of lobes in some bipolar-PNe, and demonstrates the capability of the binary-blown jets to form bipolar lobes.

The scaling of the jet velocity, $`400\mathrm{km}\mathrm{s}^1`$, is typical for jets blown by main sequence stars, which have similar escape velocities (Livio 2000). Jets blown by white dwarfs will have much higher velocities, of the order of the escape velocity from white dwarfs ($`5000\mathrm{km}\mathrm{s}^1`$). Such velocities are directly observed in some PNe and proto-PNe, e.g., a wind velocity of $`2300\mathrm{km}\mathrm{s}^1`$ was reported by Sánchez Contreras & Sahai (2001) in the proto-PNe He 3-1475. In that case, we find from equation (6) that the shocked jet material will not cool, and the bubble will be inflated as soon as the jet encounters the slow wind gas. In these cases, the lobes that are formed will have a wide opening very close to the equatorial plane, e.g., as in the symbiotic Mira He 2-104 (the Southern Crab; Corradi et al. 2001).

## 4 FAST-PROPAGATING JETS

### 4.1 Jet Interaction with the Slow Wind

When the density of the gas in the jet is much larger than that of the slow wind it will not slow much, and equation (1) is not applicable. Instead, the jet head proceeds at a velocity given by

$`v_hv_j\left[1\left({\displaystyle \frac{\rho _s}{\rho _j}}\right)^{1/2}\right];\mathrm{for}\rho _j\rho _s.`$ (12)
Most of the thermal energy released is that of the shocked slow wind material, rather than that of the shocked jet material as in the slow-propagating jet case. Thermal energy injected to the shocked slow wind material per unit time and unit length are $`\dot{E}_0=0.5A_j\rho _sv_j^3`$, and

$`e_0=\dot{E}_0/v_h0.5A_j\rho _sv_j^2,`$ (13)
respectively, where $`A_j`$ is the cross section of the jet.

As in the slow-propagating jet, bubble inflation starts when the shocked material has no time to cool. To find this distance $`r_s`$, the calculations that led to equation (6) are repeated with the following changes: ($`i`$) the energy injection rate $`0.5\dot{M}_jv_j^2`$ in equation (1) is replaced by $`\dot{E}_0`$ given above with $`A_j=4\pi \beta r^2`$, and ($`ii`$) the shocked material now is the slow wind gas. Hence the term for the total shocked mass $`\dot{M}_jt_b`$ in the derivation of equation (5) is replaced by the mass of the shocked slow wind during time $`t_b`$: $`\beta \dot{M}_st_bv_j/v_s`$, where, as before, $`t_b`$ is the expansion time of the bubble. The condition that the cooling time be longer than the expansion time $`t_b`$, becomes a condition on the location of the jet head

$`r_s6.1\times 10^{15}\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{7/3}\left({\displaystyle \frac{v_s}{15\mathrm{km}\mathrm{s}^1}}\right)^{5/6}\left({\displaystyle \frac{\dot{M}_s}{10^5M_{}\mathrm{yr}^1}}\right)^{5/6}`$
$`\times \left({\displaystyle \frac{t_b}{100\mathrm{yr}}}\right)^{1/6}\left({\displaystyle \frac{\beta }{10^3}}\right)^{1/3}\mathrm{cm}.`$ (14)
This condition is very similar to the one for slow-propagating jets (eq. 6) . As with slow-propagating jets, jets blown by main sequence stars, for which $`v_j400\mathrm{km}\mathrm{s}^1`$, will start to inflate a bubble only after propagating a significant distance from their origin; this will push slow wind material toward the equatorial plane. Jets blown by white dwarfs companions, for which $`v_j1000\mathrm{km}\mathrm{s}^1`$, on the other hand, will start to inflate a bubble very close to their origin at $`r10\mathrm{AU}`$. They will also push material toward the equatorial plane, but in a complicated flow which required 3D numerical simulations to be treated correctly (SR00).

### 4.2 The Bubble Shape

I consider two cases which bound the two extremes of a jet’s possible cross section. In the first case, the jet expands radially so that $`A_j=4\pi \beta r^2`$, at all times. In the second case, the jet expands radially at early times, but it is then recollimated by the pressure of the surrounding shocked gas or by an internal magnetic field (e.g., Kössl, Müller, & Hillebrandt 1990), such that its cross section stays constant at $`A_j=\pi a_j^2`$, where $`a_j`$ is the radius of the jet. I do not study these processes here, and simply assume that after opening up, the jet is reconfined near the place where bubble inflation starts, $`r_s10^{16}\mathrm{cm}`$ by equation (14), and I scale the jet radius by $`a_j=0.1r_s10^{15}\mathrm{cm}`$. Jets propagating with their cross section increasing at a rate much slower than $`A_jr^2`$, and even with constant cross section along a fraction of their path, are observed in different objects. Examples are the radio jet in the radio galaxy 3C 111 (Linfield & Perley 1984) and the molecular outflows in the young stellar objects HH 288 (Gueth, Schilke, & McCaughrean 2001) and HH 34 (Reipurth et al. 2000).

#### 4.2.1 Freely expanding jets

Substituting the relevant cross section, $`A_j=4\pi \beta r^2`$, in the rate of thermal energy released per unit length (eq. 13), and then in equation (7), gives the radius of the axisymmetric bubble as function of the distance, $`x`$, from the equatorial plane

$`{\displaystyle \frac{z_b}{x}}1.0\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{1/2}\left[{\displaystyle \frac{t_b(x)}{100\mathrm{yr}}}\right]^{1/2}\left({\displaystyle \frac{x}{10^{16}\mathrm{cm}}}\right)^{1/2}\left({\displaystyle \frac{\beta }{10^{10^3}}}\right)^{1/4}.`$ (15)
We find that the basic shape of the lobe is the same as in the case of a slow-propagating jet, hence, the analysis from equation (8) to the end of the previous section is applicable to the present case as well.

#### 4.2.2 Reconfined jets

Substituting the constant cross section in the thermal energy released per unit length (eq. 13), and then in equation (7), gives the radius of the axisymmetric bubble formed by a reconfined jet

$`z_b1.3\times 10^{16}\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{1/2}\left({\displaystyle \frac{a_j}{10^{15}\mathrm{cm}}}\right)^{1/2}\left[{\displaystyle \frac{t_b(x)}{100\mathrm{yr}}}\right]^{1/2}\mathrm{cm}.`$ (16)
The shape here is different from that in previous cases; rather than widening, the cylindrical lobe has a constant radius. Since the jet propagates outward, the expansion time at location $`x`$ from the equatorial plane is $`t_{\mathrm{obs}}x/v_j`$, where $`t_{\mathrm{obs}}`$ is the observing time measured from the birth of the jet . Since in the present case the jet propagate much faster than the slow-propagating case, this delay is less significant than the case analyzed in equations (10)-(11). In any case, for the same parameters as in the last equation, the bubble radius will be

$`z_b4.0\times 10^{16}\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{1/2}\left({\displaystyle \frac{a_j}{10^{15}\mathrm{cm}}}\right)^{1/2}`$
$`\times \left[\left({\displaystyle \frac{t_{\mathrm{obs}}}{1000\mathrm{yr}}}\right)0.40\left({\displaystyle \frac{x}{5\times 10^{17}\mathrm{cm}}}\right)\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^1\right]^{1/2}\mathrm{cm}.`$ (17)
The proto-PN Hen 401 has cylindrical lobes with almost constant radius of $`z_b5\times 10^{16}\mathrm{cm}`$ up to a distance of $`3.5\times 10^{17}\mathrm{cm}`$ from the central source (Sahai, Bujarrabal, & Zijlstra 1999). Farther out the lobes become narrower. Sahai et al. (1999) present detailed study of the lobes. They argue, based on the theoretical study of molecular outflows by Masson & Chernin (1993), that the structure of Hen 401 can be explained by interaction between a well collimated outflow and the slow wind, which intermediate between being radiative and being adiabatic. Some basic ingredients of their scenario apply here as well. According to the model presented here, the initial jet interaction with the slow wind is radiative, i.e., the shocked slow wind material has time to cool radiatively. This is appropriate for a jet blown by a main sequence companion; the case of a white dwarf companion is discussed below. As the jet expands to $`r10^{16}\mathrm{cm}`$ (eq. 14 above) the adiabatic phase starts, with the inflation of a bubble (which will become the observed lobes; one lobe at each side of the equatorial plane). This part of the bubble compresses the material near the equatorial plane; this explains the wide opening of the lobe close to the equatorial plane, which cannot be explained with a fully radiative flow. The jets then expand with a constant cross section. Further out, the jets may become weaker, and the interaction becomes unstable (Sahai et al. 1999). To quantitatively match the structure of Hen 401, $`z_b=5\times 10^{16}\mathrm{cm}`$ up to $`x=3.5\times 10^{17}\mathrm{cm}`$, the age of the jet and its velocity should be known. This information is not available, hence, I consider a few examples. We can take $`a_j=3\times 10^{15}\mathrm{cm}`$, keeping the other parameters as in equation (17); or we can take $`t_{\mathrm{obs}}=1500\mathrm{yr}`$, $`v_j=500\mathrm{km}\mathrm{s}^1`$, and keeping $`a_j`$ as in equation (17). For a white dwarf companion, we set $`v_j5000\mathrm{km}\mathrm{s}^1`$ and take $`a_j2\times 10^{14}\mathrm{cm}13\mathrm{AU}`$, keeping the observing time as in equation (17). In this case, the bubble inflation starts much closer to the binary system (eq. 14). To summarize, the flow structure presented here, despite its several assumptions, e.g., neglecting the decrease in the slow wind density along the $`z`$ direction (perpendicular to the symmetry axis), can explain the structure of the lobes in Hen 401, and similar PNe having cylindrical lobes (see Sahai et al. 1999 for more examples).

## 5 SUMMARY

In the present paper, I studied the interaction of jets, or a collimated fast wind (CFW), with a previously blown AGB wind, and reproduced some basic properties of bipolar PN lobes. In reality, the jets (or CFW) are blown by an accreting compact companion simultaneously with the slow wind. To facilitate analytical solution, I had to assume axisymmetric flow, as well as several other simplifying assumptions. In addition, I studied only two extreme cases, a slow-propagating and fast-propagating jets, which occur when the slow wind density is much larger, or much lower, than the jet density, respectively. The goal was to derive simple expressions from which the properties of PN lobes blown by jets can be inferred. This can be used both to analyze observations (see discussions following eqs. 6 and 17), and to guide future numerical simulations.

The main results can be summarized as follows:
(1) The interaction of jets (or CFW) blown by main sequence stars, having a speed of $`v_j300500\mathrm{km}\mathrm{s}^1`$, with the slow AGB stellar wind close to the central binary system is radiative, i.e., the shocked gas cools on a short time. Further away, the interaction becomes adiabatic, i.e., a bubble of hot gas which cools slowly is formed. The distance from the center where this transition takes place for the two types of jets studied here is given by equations (6) and (14). The expressions are scaled for jets blown by main sequence companions. Jets blown by white dwarf companions move much faster, $`v_j5000\mathrm{km}\mathrm{s}^1`$, and start to inflate a bubble very close to the binary system.
(2) As noted by SR00, the interaction between the CFW and the slow wind pushes material toward the equatorial plane. The formation of the inflated bubble at a distance from the equatorial plane in jets blown by main sequence stars, increases the efficiency of this process.
(3) The significance of the transition from the radiative to the adiabatic phase to the structure of the lobes, implies that radiative cooling must be incorporated in numerical codes intended to study the formation of bipolar PNe.
(4) The basic interaction of a radially expanding jet with the slow wind forms an axisymmetric bubble which is narrow close to the central star(s) and far from the central star(s); its widest diameter is somewhere in-between. This type of structure is observed in many bipolar PNe.
(5) If a jet is recollimated, e.g., by magnetic fields inside the jet, such that its cross section stays constant, then the lobe widens close to the center, but then acquires a cylindrical shape with a constant diameter. Such lobes are observed in several bipolar PNe.

I thank Liz Blanton for comments on the original manuscript. This research was supported in part by grants from the US-Israel Binational Science Foundation.