text:
# Comparison of source detection procedures for XMM-Newton images

## 1 Introduction

X-ray astronomy has entered a new era now that Chandra and XMM-Newton are in orbit. Their high sensitivities and unprecedented image qualities bear great promises but also pose new challenges. In this paper, we outline problems of object detection in X-ray images that were not previously encountered. In doing so, we compare the performances of various detection techniques on simulated XMM-Newton test images, incorporating the main instrumental characteristics.

The X-ray observations consist of counting incoming photons one by one, recording their time of arrival, position and energy. Later, the event list is used to create images for a given pixel scale and energy band. Various X-ray telescope effects complicate this simple picture – the point spread function (PSF) and the telescope effective area (the vignetting effect), both dependent on the off-axis angle and incoming photon energy; detector effects like quantum efficiency variations, different zones not exposed to X-ray photons; environmental and background effects like solar flares and particle background. Even for relatively large exposures, the X-ray images could contain very few photons, and some sources could contain only a few tens of photons spread over a large area. Consequently, it is important for the source detection and characterization procedures to be able to cope with these difficulties.

As an example, the same hypothetical input situation is shown schematically in Fig. 1 for ROSAT 111http://wave.xray.mpe.mpg.de/rosat, XMM-Newton 222http://xmm.vilspa.esa.es/ and Chandra 333http://chandra.harvard.edu/ . XMM-Newton ’s rather large PSF, coupled with its higher sensitivity, leads to the detection of more objects but also to blending and source confusion, which become severe for long exposures depending on the energy band. Confusion problems in the hard band above 5 keV are less important, given the smaller number of objects and smaller count rate of energetic photons. Thus, why we concentrate our analysis mainly on source detection problems for the more complicated case of the XMM-Newton energy bands below 5 keV.

Each X-ray mission provides data analysis packages – EXSAS for ROSAT (Zimmermann et al. exsas (1994)), CIAO for Chandra (Dobrzycki et al. ciao (1999)) and XMM-Newton Science Analysis System (XMM-SAS444http://xmm.vilspa.esa.es/sas/). They include procedures for source detection, and in this paper we estimate and compare their performances on simulated images using various types of objects. These procedures make use of techniques such as Maximum Likelihood (ML), Wavelet Transformation (WT), Voronoi Tessellation and Percolation (VTP).

In Section 2 we describe the X-ray image simulations. A short presentation of the detection procedures is given in Sec. 3. Tests using only point sources are presented in Sec. 4, and extended sources in Sec. 5. We have analyzed realistic simulations of a shallow and a deep extragalactic field with only point sources in Sec. 6 and with extended objects in Sec. 7 for an exposure of 10 ks. Finally, we investigate the problems of confusion and completeness in two energy bands – \[0.5-2\] and \[2-10\] keV for two exposures – 10 ks and 100 ks (Sec. 8). Sec. 9 presents the conclusions. (H$`{}_{0}{}^{}=50`$ km/s/Mpc, $`h=0.5`$, q$`{}_{0}{}^{}=0.5`$ and $`\mathrm{\Omega }_0=1`$ are used).

## 2 Simulation of X-ray images

The simulations are essential to understand and qualify the behavior of the different detection and analysis packages. We have developed a simulation program that generates X-ray images for given exposure times with extended and point-like objects. It takes into account the main instrumental characteristics of XMM-Newton and the total sensitivity of the three EPIC instruments. The procedure is fast and flexible and is made of two independent simulation tasks: object generation (positions, fluxes, properties) and instrumental effects. A possibility to apply the instrumental response directly over images is also implemented, especially useful when one wants to use sky predictions from numerical simulations (cf. Pierre et al. mp00 (2000)).

A summary of the simulated images parameters are given in Tab. 1.

The point-like sources are assumed to be AGNs or QSOs with a power law spectrum with a photon index of 2 and flux distribution following the $`\mathrm{log}N\mathrm{log}S`$ relations of Hasinger et al. (has98 (1998), has01 (2001)) and Giacconi et al. (gia01 (2000)) in the two energy bands.

The PSF model is derived from the current available calibration files555http://xmm.vilspa.esa.es/ccf/. On-board PSF data is generally in very good agreement with the previous ground based calibrations (Aschenbach et al. asch00 (2000)). We must stress that the model PSF is an azimuthal average and in reality, especially for large off-axis angles, its shape can be very distorted. However, in the analytic model (Erd et al. calbook (2000)), the off-axis and energy dependences are not available yet. This is not crucial, as the energy dependence in the bands used is moderate and we confine all the analysis inside $`10\mathrm{}`$ from the centre of the field-of-view where the PSF blurring is negligible.

The extended objects are modeled by a $`\beta `$-profile (Cavaliere & Fusco-Femiano cff76 (1976)) with fixed core radius $`r_c=250h^1`$ kpc and $`\beta =0.75`$. A thermal plasma spectrum (Raymond & Smith rs77 (1977)) is assumed for different temperatures, luminosities and redshifts.

The source spectra (extended and point-like) are folded with the spectral response function for the total sensitivity of the three XMM-Newton EPIC instruments (MOS1, MOS2 and pn with thin filters) by means of XSPEC (Arnaud xspec (1996)) to produce the count rates in different energy bands. The actual choice of the energy bands is not important for this comparison study, although some objects can be more efficiently detected in particular energy ranges.

As an example, we show in Tab. 2 the resulting count rates for extended sources assuming that they represent an average cluster of galaxies.

The background in the simulations includes realistic cosmic diffuse and background values derived from the XMM-Newton in-orbit measurements in the Lockman Hole (Watson et al. wat00 (2001)).

The calculated count rates for the objects and the photons of the background are subject to the vignetting effect – some photons are lost due to the smaller telescope effective area at given off-axis angle, depending on the incoming photon’s energy. We have parametrized the vignetting factor – the probability that a photon at an off-axis angle $`\theta `$ to be observed – as polynomials of fourth order in two energy bands: \[0.5-2\] and \[2-10\] keV, using the latest XMM-Newton on-flight calibration data. For example, a photon at $`\theta =10\mathrm{}`$ has a 53% chance of being observed in \[0.5-2\] keV and 48% in \[2-10\] keV.

Further instrumental effects such as quantum efficiency difference between the CCD chips, the gaps between the chips, out-of-time events, variable background, pile-up of the bright sources are not taken into account – their inclusion is not relevant for our main objective.

## 3 Detection procedures

Without attempting to provide a review of the available techniques in the literature, we briefly describe here the procedures we have tested. They are summarized in Tab. 3.

### 3.1 Sliding cell detection and maximum likelihood (ML) method

Historically, the sliding cell detection method was first used for Einstein Observatory observations (e.g. EMSS – Gioia et al. emss (1990)). It is included in ROSAT , Chandra and XMM-Newton data analysis tools and a good description can be found in the specific documentation for each of those missions.

The X-ray image is scanned by a sliding square box and if the signal-to-noise of the source centered in the box is greater than the specified threshold value it is marked as an object. The signal is derived from the pixel values inside the cell and noise is estimated from the neighboring pixels. Secondly, the objects and some zone around them are removed from the image forming the so-called “cheese” image which is interpolated later by a suitable function (generally a spline) to create a smooth background image. The original image is scanned again but this time using a threshold from the estimated background inside the running cell to give the map detection object list.

The procedure is fast and robust and does not rely on a priori assumptions. However it has difficulties, especially in detecting extended features, close objects and sources near the detection limit. Many refinements are now implemented improving the sliding cell method: (1) consecutive runs with increasing cell size, (2) matched filter detection cell where the cell size depends on the off-axis angle. However, the most important improvement was the addition of the maximum likelihood (ML) technique to further analyze the detected sources.

The ML technique was first applied to analyze ROSAT observations (Cruddace et al. ml88 (1988), cru91 (1991), Hasinger et al. has93 (1993)). It was used to produce all general X-ray surveys from ROSAT mission (e.g. RASS – Voges et al. rass (1999), WARPS survey – Ebeling et al. warps (2000)). The two lists from local and map detection passes can be merged to form the input objects list for the ML pass. It is useful to feed the ML procedure with as many candidate objects as possible, having in mind that large numbers of objects could be very CPU-expensive. The spatial distribution of an input source is compared to the PSF model – the likelihood that both distributions are the same – is calculated with varying the input source parameters (position, extent, counts) and the corresponding confidence limits can be naturally computed. A multi-PSF fit is also implemented which helps in deblending and reconstructing the parameters of close sources. In the output list, only sources with a likelihood above a threshold are kept.

The ML method performs well and has many valuable features, however, it has some drawbacks – it needs a PSF model to perform the likelihood calculation and thus favours point-like source analysis, extent likelihood could be reliably taken only for bright sources, it cannot detect objects which are not already presented in the input list (e.g. missing detections in the local or map passes).

Here we have used EMLDETECT – an implementation of the method specifically adapted for XMM-SAS (Brunner b96 (1996)). In the map mode sliding cell pass we used a low signal-to-noise ratio ($`3\sigma `$) above the background in order to have as many as possible input objects for the ML pass. The likelihood limit (given by $`L=\mathrm{ln}P`$, where $`P`$ is the probability of finding an excess above the background) was taken to be 10, which corresponds roughly to $`4\sigma `$ detection. A multi-PSF fitting mode with the maximum of 6 simultaneous PSF profile fits was used.

### 3.2 VTP

VTP – the Voronoi Tessellation and Percolation method (Ebeling & Wiedenmann vtp1 (1993), Ebeling vtp2 (1993)) is a general method for detecting structures in a distribution of points (photons in our case) by choosing regions with enhanced surface density with respect to an underlying distribution (Poissonian in X-ray images). It treats the raw photon distribution directly without any recourse to a PSF model or a geometrical shape of the objects it finds. Each photon defines a centre of a polygon in the Voronoi tessellation image and the surface brightness is simply the inverse area of the polygon (assuming one single photon per cell). The distribution function of the inverse areas of all photons is compared to that expected from a Poisson distribution and all the cells above a given threshold are flagged and percolated, i.e. connected to form an object. This method was successfully used with ROSAT data (Scharf et al. s97 (1997)) and is currently incorporated in the Chandra DETECT package (Dobrzycki et al. ciao (1999)).

Apart from these advantages, VTP has some drawbacks which are especially important for XMM-Newton observations: (1) because of the telescope’s high sensitivity and rather large PSF with strong tails, the percolation procedure tends to link nearby objects; (2) excessive CPU time for images with relatively large number of photons; (3) there is no simple way to estimate the extension of objects.

### 3.3 Wavelet detection procedures

In the past few years a new approach has been extensively used: the wavelet technique (WT). This method consists in convolving an image with a wavelet function:

$$w_a(x,y)=I(x,y)\psi (\frac{x}{a},\frac{y}{a}),$$
(1)
where $`w_a`$ are the wavelet coefficients corresponding to a scale $`a`$, $`I(x,y)`$ is the input image and $`\psi `$ is the wavelet function. The wavelet function must have zero normalization and satisfy a simple scaling relation

$$\psi _a(x,y)=\frac{1}{a}\psi _1(\frac{x}{a},\frac{y}{a}).$$
(2)
The choice of $`\psi `$ is dictated by the nature of the problem but most often the second derivative of the Gaussian function or the so called “Mexican hat” is used.

The WT procedure consists of decomposing the original image into a given number of wavelet coefficient images, $`w_a`$, within the chosen set of scales $`a`$. In each wavelet image, features with characteristic sizes close to the scale are magnified and the problem is to mark the significant ones, i.e. those which are not due to noise. In most cases, this selection of significant wavelet coefficients cannot be performed analytically because of the redundancy of the WT introducing cross-correlation between pixels. For Gaussian white noise, $`w_a`$ are distributed normally, allowing easy thresholding. This is not the case for X-ray images which are in the Poissonian photon noise regime.

Various techniques were developed for selecting the significant wavelet coefficients in X-ray images. In Vikhlinin et al. (vik97 (1997)) a local Gaussian noise was assumed; Slezak et al. (sle94 (1994)) used the Ascombe transformation to transform an image with Poissonian noise into an image with Gaussian noise; in Slezak et al. (sle93 (1993)), Starck & Pierre (sp98 (1998)) a histogram of the wavelet function is used. In recent years a technique based on Monte Carlo simulations is used successfully (e.g. Grebenev et al. gre95 (1995), Damiani et al. dam97 (1997), Lazzati et al. l99 (1999)).

Once the significant coefficients at each scale are chosen, the local maxima at all scales are collected and cross-identified to define objects. Different characteristics, such as centroids, light distribution etc., can be computed, as well as an indication of the source size at the scale where the object wavelet coefficient is maximal.

WT has many advantages – the multiresolution approach is well suited both for point-like and extended sources but favours circularly symmetric ones. Because of the properties of the wavelet function a smoothly varying background is automatically removed. Extensive description of wavelet transform and its different applications can be found in Starck et al. (sta98 (1998)).

In this work we have tested two WT procedures:

one of the Chandra wavelet-based detection techniques (Freeman et al. fre96 (1996), Dobrzycki et al. ciao (1999)). It uses a “Mexican hat” as wavelet function and identifies the significant coefficients by Monte Carlo simulations. The background level needed to empirically estimate the significance is taken directly from the negative annulus of the wavelet function. Source properties are computed inside the detection cell defined by minimizing the function $`|\mathrm{log}r_{PSF}\mathrm{log}\sigma _F|`$, where $`r_{PSF}`$ is the size of the PSF encircling a given fraction of the total PSF flux and $`\sigma _F`$ is the size of the object at the scale closest to the PSF size. No detailed PSF shape information is needed to perform this minimization – just $`r_{PSF}`$ as a function of the off-axis angle. The control parameters which need attention are the significance threshold and the set of scales. We have used a significance threshold of $`10^4`$ corresponding to 1 false event in 10000 ($`4\sigma `$ in Gaussian case) and “$`\sqrt{2}`$ sequence” for the scales, where $`a=1,\sqrt{2},2,2\sqrt{2},4,3\sqrt{2}\mathrm{}16`$.

XMM-SAS package, based on a wavelet analysis. It is very similar to WAVDETECT but implements some new ideas. The sources are assumed to have a Gaussian shape in order to analytically derive their counts and extent. Currently, the PSF information is ignored, which can be regarded as a serious drawback. We have used the significance threshold of $`10^4`$ ($`4\sigma `$) and wavelet scales 1, 2, 4, 8 and 16.

### 3.4 Mixed approach

This method combines a source detection by an elaborated procedure over filtered/smoothed raw photon images.

The use of such a mixed approach is motivated by the fact that procedures for source detection in astronomical images have been developed for many years and the steps and problems of deblending, photometry, classification of objects are now quite well understood. The raw photon image manipulations can be performed with very simple smoothing procedures (for example a Gaussian convolution) or with more sophisticated methods like the “matching filter” technique, adaptive smoothing or multiresolution (wavelet) filtering.

We have used two different types of raw image filtering:

Gaussian convolution. For our simulated images we applied two convolutions with FWHM=$`12\mathrm{}`$ and $`20\mathrm{}`$ constrained by the image characteristics (see Sec. 2).

Multiresolution iterative threshold filtering (MR/1  package, Starck et al. sta98 (1998)) by means of “à trous” (with holes) wavelet method and Poissonian noise model (Slezak et al. sle93 (1993), Starck & Pierre sp98 (1998)) to flag the significant wavelet coefficients (also known as auto-convolution or wavelet function histogram method). The control parameters are the significance threshold, wavelet scales and the number of iterations for the image reconstruction. We took $`10^4`$ significance threshold (corresponding to $`4\sigma `$ for Gaussian distribution), 6 scales for the “à trous” method allowing analysis of structures with characteristic sizes up to $`2^6`$ pixels ($`4\mathrm{}`$). We have used 25 iterations or a reconstruction error less than $`10^5`$. More iterations can improve the photometry but can also lead to the appearance of very faint real or sometimes false features. Optionally the last-scale wavelet smoothed image in the iterative restoration process is used for analysis of large-scale background variations or very extended features.

Source detection over the filtered images was performed with SExtractor (Bertin & Arnouts sex (1996)) – one of the most widely used procedures for object detection in astronomy with a simple interface and very fast execution time. Originally, SExtractor was developed for optical images but it can be applied very successfully on X-ray images after they are properly filtered. The flexibility of SExtractor is ensured by a large number of control parameters requiring dedicated adjustments. In the following we will briefly outline those having a significant influence on our results:

* Detection and analysis threshold – important parameters for Gaussian convolved images. It is preferable to use a value greater than $`4\sigma `$, even sacrificing very faint sources but reducing the number of false detections. For MR/1  filtering this is irrelevant because the features in the images are already at $`4\sigma `$ significance.
* Background mesh size – influences object detection as well as photometry. Using too coarse a background grid (less grid points) will lead to a smoother background image but some faint objects will be missed in the detection. While too detailed a map (more grid points) will lead to many false detections and worse photometry, depending on the local noise properties. We have found after many experiments that local background meshes with 32 to 64 nodes give best results balancing both effects.
* Minimum detection area – in some cases this parameter can help to avoid spurious detections especially when one uses wavelet scales below the characteristic size of objects.
* Deblending parameters – the number of the deblending levels and the minimum contrast. We used the extreme values of 64 levels and zero contrast so that any saddle point will lead to object splitting.
* Photometry – we have decided to use the automatic aperture photometry procedure implemented in SExtractor following Kron (kro80 (1980)) and Infante (inf87 (1987)) which gives better results than the fixed aperture or isophotal methods.
* Classification – we discuss the classification in the next subsection.

### 3.5 Classification

The problem of classifying sources as resolved or unresolved (or stars/galaxies, extended/point-like) has a long history and discussions can be found in widespread detection packages like FOCAS (Valdes et al. focas (1990)), INVENTORY (Kruszewski inventory (1988)), SExtractor (Bertin & Arnouts sex (1996)), Neural Extractor (Andreon et al. next (2000)). This task becomes even more difficult for X-ray images – we have already mentioned various X-ray telescope effects that blur and distort the images as a function of the off-axis angle.

– the extension likelihood is calculated by the deviation of the object fitted from the PSF model when varying the size (extension) of the object. It was used with ROSAT observations to flag possible extended objects and in most of the cases the results were positive. The extension likelihood depends on the chosen PSF model and for faint objects is not reliable (De Grandi et al. sdg97 (1997)).

– the ratio between object size and the PSF size ($`R_{PSF}`$) is used for classification purposes. It depends on the adopted definition for object size and can be misleading, especially for faint objects.

– one parameter for classification is the stellarity index. It is based on neural-network training performed over set of some $`10^6`$ simulated optical images. There are 10 parameters included in the neural network training – 8 isophotal areas, the maximal intensity and the “seeing”. The only controllable parameter is the “seeing” which, however, has no meaning for space observations. It can be tuned to the mean PSF size.
The half-light radius $`R_{50}`$ is another classification parameter. It can be used as a robust discriminator because it depends only on the photon distribution inside the object (“luminosity profile”), which is assumed to be different for resolved and unresolved objects.

Additional indications might be taken into account when the classification is ambiguous – e.g. the wavelet scale at which object’s wavelet coefficient is maximal or even spectral signatures.

## 4 Test 1 – point-like sources

### 4.1 Input configuration

We address the problem of point-like sources separated by $`15\mathrm{}`$ (half-energy width of the on-axis PSF), $`30\mathrm{}`$ and 60$`\mathrm{}`$ with different flux ratios. We include the PSF model and background but do not apply the vignetting effect.

The raw input test image is shown in Fig. 2 together with its Gaussian convolution, MR/1  wavelet filtering and WAVDETECT output image. Visually, the Gaussian image is quite noisy, while there are few spurious detections in the WT images.

### 4.2 Detection rate and positional errors

The number of missing detections and false objects are shown in Tab. 4.

The one sigma input-detect position differences are less than the FWHM of the PSF ($`6\mathrm{}`$) for all procedures and the maximum occurs for the blended objects, as expected. Note the large number of spurious detections with WAVDETECT, VTPDETECT and EMLDETECT.

### 4.3 Photometric accuracy

The results for the photometry in terms of the inferred to the input counts are shown in Fig. 3.

Only EMLDETECT detects two of the six fainter objects. None of the other procedures separates the objects and consequently the inferred counts are a blend from both sources.

The proximity of objects influences the detection and the photometry. EMLDETECT and EWAVELET show the best detection rate results while all the other procedures miss one of the faintest objects.

We can safely assume that the objects are well separated. The recovery of the properties is informative for the performance of the tested procedures. It is clear that the general flux reconstruction error (taken to be the spread of the points around the unity line in Fig. 3) is about 15% for brighter sources and goes down to 20-25% for the faintest ones. \[In our 10 ks exposure tests and with the adopted background in band \[0.5-2\] keV, we assume the objects with input counts of 20 photons ($`10^{15}`$ erg/s/cm2) to be at the detection limit when there is no confusion by nearby sources.\]

### 4.4 Discussion

After this simple test we can eliminate the VTPDETECT: in addition to the very large execution time, some of the VTP-detected object centres were shifted by more than $`20\mathrm{}`$ from their input positions – a consequence of its ability to detect sources with different shapes where the object center can be far from the input position. Moreover, VTPDETECT percolates all the double sources into single objects at $`\delta r=30\mathrm{}`$, which all other procedures were able to separate.

No procedure unambiguously shows best results – both in terms of the detection rate, spurious sources and photometric reconstruction. EMLDETECT outperforms the others in terms of detection rate but with the price of many spurious detections. Using exactly the same PSF model as the one hard-coded in EMLDETECT leads to much better photometric reconstruction.

All other procedures are comparable: EWAVELET showing better detection but its photometric reconstruction is far from satisfactory – about half of the photons were lost at $`\delta r=30\mathrm{}`$ and $`60\mathrm{}`$, because of the assumed Gaussian shape used to derive analytically extension and counts. We have applied a simple correction for the object size to arrive at the good photometric results for EWAVELET presented on Fig. 3.

## 5  Test 2 – point-like plus extended objects

### 5.1 Input configuration

This test is similar to Test 1, but we have replaced some of the point-like sources by extended ones generated as described in Sec. 2. The raw photon image with input counts indicated and its representations are shown in Fig. 4.

### 5.2 Detection rate and positional errors

The number of missed and false detections are shown in Tab. 5. An increase of the searching radius to $`20\mathrm{}`$ was needed: at $`\delta r=15\mathrm{}`$ the blending tends to shift the centroid towards the point-like source. Note that this situation is a clear case for source confusion: if we take the closest neighbour (the point-like source in some cases) as the cross-identification from the input list, we shall overestimate the flux more than two-fold, while the true representation is the extended object.

Some changes were needed for the procedures not based on the wavelet technique in order to avoid splitting of the bright extended objects into sub-objects: increase of the Gaussian convolution FWHM to $`20\mathrm{}`$, and multi-PSF fit for EMLDETECT. In the Gaussian case, the larger smoothing length smears some of the point-like sources, leading to non-detection. EMLDETECT splitting persists even with the maximum number of the PSFs fitted to the photon distribution (currently it is capable of simultaneously fitting up to 6 PSFs).

### 5.3 Photometric reconstruction

The inferred-to-input source counts ratio $`SCTS(out)/SCTS(in)`$ is shown in Fig. 5. We will not consider EWAVELET as its photometric results and detection of extended sources were quite unsatisfactory. The simple correction technique (as in Test 1) based on the PSF does not work – the extended objects profile has different shape from the PSF.

Only MR/1+SE gives good results for the flux restoration of the extended objects, but overestimating the counts for the faintest one. WAVDETECT misses almost half of the input photons while EMLDETECT splitting leads to very poor results.

All procedures give bad restoration results with MR/1+SE performing best again for the extended sources. The proximity of the objects leads to an overestimation of the point-like source counts and an underestimation of the extended object counts. There is no simple way to correct for this effect, but can be done using a rather elaborated iterative procedure involving extended object profile fitting.

Point-like source results are relatively similar with all procedures – the source counts are slightly overestimated due to the extended object halo even at $`60\mathrm{}`$. The problems of WAVDETECT and EMLDETECT and the recovery of the extended objects counts are quite obvious. Again MR/1+SE is the best performing procedure with extended objects flux uncertainty about 25-30%.

### 5.4 Object classification

An important test is the ability of the procedures to classify objects and to allow further analysis of complicated cases of blending. The MR/1+SEresults for the clasification by means of the half-light radius and stellarity index (c.f. Sec. 3.5) are shown in Fig. 6, overlayed over the results from 10 simulated images with only point-like sources. Clearly the detected extended objects with MR/1+SE fall into zones not occupied by point-like sources. Note, however, that the two detected point-like sources at $`30\mathrm{}`$ will be mis-clasified as extended objects – the proximity not only influences the source counts, overestimated by more than 4-5 times (Fig. 5), but also the object profile and consequently the classification.

Fig. 7 shows the WAVDETECTclassification – the ratio of the object size to the PSF size ($`R_{PSF}`$). The results are more ambiguous with WAVDETECT (Fig. 7) compared to MR/1+SE. The results with EMLDETECT and its classification parameter (extension likelihood) were very unsatisfactory due to the extended object splitting. More comprehensive discussion of the simulations and the classification is left for Sec. 7.

### 5.5 Discussion

Clearly EMLDETECT and WAVDETECT have problems in restoring the fluxes of extended objects. We have already discussed the splitting difficulties of EMLDETECT. The explanation for WAVDETECT’s bad results is that the wavelet scale at which the detected object size is closer to the PSF size defines the source detection cell (in which the flux is computed). If the characteristic size of an object is larger than the PSF size (i.e. an extended object) this procedure will tend to underestimate the flux.

We can safely accept the MR/1+SE procedure as the best performing for detection and characterization both for point-like and extended objects. We must stress however that one cannot rely on the flux measurements when there are extended and point-like sources separated by less than $`30\mathrm{}`$. The proximity affects also the classification of the point-like sources. Using the classification and then performing more complicated analysis like profile fitting and integration for the extended sources can improve a lot the restoration. In realistic situations we can expect very often problems of this kind, especially with XMM-Newton .

## 6 Test 3 – Realistic model with only point-like sources

### 6.1 Input configuration

We simulate an extragalactic field including only point-like sources with fluxes drawn from the $`\mathrm{log}N\mathrm{log}S`$ relation (Hasinger et al. has98 (1998), has01 (2001), Giacconi et al. gia01 (2000)). PSF, vignetting and background models are applied as described in Sec. 2. The aim is to test the detection procedures in more realistic cases where confusion and blending effects are important and not controlled. The raw photon image is shown in Fig. 8 together with its visual representation – the same input configuration for a much larger exposure time and no background, only keeping the objects with counts greater than 10. It displays better the input object sizes, fluxes and positions and it is instructive to compare it to the MR/1  filtered and WAVDETECT images shown on the same figure.

### 6.2 Cross-identification and positional error

We need to define a searching radius in order to cross-identify the output and the input lists. The input list contains many objects with counts well below the detection limit ($`\mathrm{log}N\mathrm{log}S`$ extends to very faint fluxes) and a lower limit must be chosen. For each detected object, we search for the nearest neighbour inside a circle within the reduced input list.

The positional difference for the brightest detected sources (more than 100 counts) in the inner $`10\mathrm{}`$ from the center of the FOV is shown in Tab. 6. The region beyond $`10\mathrm{}`$ is subject to serious problems caused by the vignetting and PSF blurring, the detected object centroid can be few PSFs widths from the true input identification.

We therefore adopt the following cross-identification parameters: the input list is constrained to counts greater than 10; a $`6\mathrm{}`$ searching radius; we consider only the central $`10\mathrm{}`$ of the FOV.

### 6.3 Detection rate and photometric reconstruction

The detection rate and flux reconstruction results are shown in Fig. 9. There are different effects playing a role in the distribution and the numbers of missed and false detections:

“false” detections – non-existent objects, or two or more sources blended into a single detected object. The result will be a “false” detection if the blended objecs are not in the input list (count limit) or the merged object centroid is beyond the searching radius.

source confusion – in the cross-identification process the nearest neighbour to the detected source is not the true assignment; or as in case (1), when a blend of objects is wrongly identified by one input source.

missed detections – depending on the local noise properties, some objects can be missed even if their input counts are above the adopted limiting counts for cross-identification.

### 6.4 Discussion

The results in terms of detection rate are similar for all procedures. The best detection rate shows G+SE  but at the price of twice as many false detections.

The photometry reconstruction for the sources above 50 counts shows a spread about 10-15% for the WT based methods and $`25\%`$ for EMLDETECT. However, EMLDETECT clearly outperforms the other procedures when we use the same PSF model as the one hard-coded into the programme. This fact shows that using a correct PSF representation has a crucial importance for the ML technique. More discussion about the detection limits, completeness and confusion is left for Sec. 8.

## 7 Test 4 – sky models with point-like and extended objects

### 7.1 Input configuration

We have chosen exactly the same input point-like sources configuration as in Test 3 and we have added 5 extended objects. They may be regarded as clusters of galaxies at different redshifts with the same $`\beta `$-profiles and moderate X-ray luminosity $`L_{[210]keV}3\times 10^{44}`$ erg/s (c.f. Tab. 2). Now the objective is to estimate the capabilities of the procedures to detect and identify extended objects. The difference from Test 2 is the arbitrary positions of the point-like sources leading to different local and global background properties and to uncontrolled confusion effects.

The input configuration and wavelet filtered and output images are given in Fig. 10.

### 7.2 Positional and photometric reconstruction

Results for the point-like sources will not be presented, because the input configuration (position, $`\mathrm{log}N\mathrm{log}S`$, background) is exactly the same as in the previous test. It was shown in Test 2 that the presence of a point source even with moderate counts in the vicinity of a faint extended source could lead to confusion and even non-detection. Of course, the presence of a cluster will affect the detection and photometry of the point-like objects in its vicinity, but this effect is of minor concern for this test and has been already discussed (Test 2).

The detection rate, input-detect position offsets, detected counts and detected-to-input counts ratio are shown in Tab. 7.

As to the positional errors, it was already concluded that the centres of the extended object can be displaced by more than the adopted searching radius for point-like sources (Sec. 5). The differences in positions shown in Tab. 7, especially for fainter objects, are 3-4 times larger than the one sigma limit for point-sources inside the inner $`10\mathrm{}`$ of the FOV (Tab. 6).

It is confirmed again that EMLDETECT and WAVDETECT are not quite successful in charactering extended objects. But note that this time the results for MR/1+SE and G+SE  are worse than the results in Test 2 – the rate of lost photons being about 20-30%. Also, the flux of the clusters at $`z=1.5`$ and $`2`$ is overestimated, suggesting blending with faint nearby point-like sources.

### 7.3 Classification

To classify the detected objects we have performed many simulations with only point-like sources as in Test 3. Ten simulated images were generated with the same $`\mathrm{log}N\mathrm{log}S`$ and background, but with different and arbitrary positions of the input sources. Exactly the same parameters were used for filtering, detection and characterization.

Two classification parameters were used: the half-light radius ($`R_{50}`$) and the stellarity index from SExtractor based procedures. The results are shown in Fig. 11. We do not show results with WAVDETECT and EMLDETECT classification parameters because their unsatisfactory results were confirmed, as in Test 2.

We can see the excellent classification based on the stellarity index and half-light radius: in the inner $`10\mathrm{}`$, for objects with more than 20 detected counts, stellarity less than 0.1 and $`R_{50}20\mathrm{}`$ we have 15 incorrect assignments from 1287 detections ($`1\%`$).

## 8 Test 5 – completeness and confusion

In this section we investigate the confusion and completeness problems for XMM-Newton shallow and deep observations like the first XMM Deep X-ray Survey in the Lockman Hole (Hasinger et al. has01 (2001)).

A set of 10 images with exposure times of 10 ks and 100 ks in the energy bands \[0.5-2\] and \[2-10\] keV were generated; the fluxes were drawn using the latest $`\mathrm{log}N\mathrm{log}S`$ relations from Hasinger et al. (has01 (2001)) and Giacconi et al. (gia01 (2000)). Detection and analysis were performed with exactly the same parameters for all simulations: detection threshold, analysis threshold, background map size, detection likelihood, etc. (see Sec. 3). Cross-identification was achieved using the input sources above 10 counts and 30 counts for 10 ks and 100 ks exposures respectively. Lowering the count limits yields more cross-identifications but increases considerably the number of spurious detections.

The input image for 100 ks and \[0.5-2\] keV band is shown in Fig. 12. The inner $`10\mathrm{}`$ zone where all analysis is performed is indicated, as well as the total XMM-Newton field-of-view. It is informative to compare it with the images for 10 ks in Fig. 8.

In order to estimate the effect of confusion we have generated images with only point-like sources, distributed on a grid such to avoid close neighbours, and with fixed fluxes spanning the interval $`[10^{16},10^{13}]`$ erg/s/cm2.

In the following we discuss some important points.

* Confusion and completeness
The detection rate of the input sources as a function of flux (Fig. 13) indicates that confusion problems are significant for 100 ks exposures in the \[0.5-2\] keV band. They are marginal for 100 ks exposures in \[2-10\] keV band and absent for 10 ks in both bands.

Energy conversion factors and the limiting fluxes below which more than 10% of the input objects are lost are shown in Tab. 8.
* Differential flux distribution
The differential flux distributions for 100 ks in \[0.5-2\] keV obtained by MR/1+SE and EMLDETECT are shown on Fig. 14. Spurious detections appears to be numerous with EMLDETECT below $`5\times 10^{16}`$ erg/s/cm2 and tend to compensate the loss of sensitivity and confusion. MR/1+SE appears to be less affected and allows us to set a conservative flux limit of $`6\times 10^{16}`$ erg/s/cm2($`90`$ photons for 100 ks), below which the incompleteness becomes important – 65% of the input sources are lost between $`3\times 10^{16}`$ and $`6\times 10^{16}`$ erg/s/cm2.
* $`\mathrm{log}N\mathrm{log}S`$
The $`\mathrm{log}N\mathrm{log}S`$ functions in \[0.5-2\] keV for 10 ks and 100 ks exposures are shown in Fig. 15. The inferred $`\mathrm{log}N\mathrm{log}S`$ (by simple source counting) are in very good agreement with the input ones down to fluxes about $`2\times 10^{15}`$ and $`6\times 10^{16}`$ erg/s/cm2 with MR/1+SEand $`10^{15}`$ and $`3\times 10^{16}`$ erg/s/cm2 for EMLDETECTfor 10 ks and 100 ks respectively. Although there are confusion and completeness problems for 100 ks starting at $`56\times 10^{16}`$ erg/s/cm2as discussed above, their effect is completely masked in the EMLDETECT integral distribution, which seems to be in very good agreement down to $`3\times 10^{16}`$ erg/s/cm2(the lower flux limit for the Lockman Hole Deep Survey analysis of Hasinger et al. has01 (2001)).
* Photometric accuracy
The photometric reconstruction is relatively similar for 10ks and 100ks in for fluxes greater than $`2\times 10^{15}`$ and $`6\times 10^{16}`$ erg/s/cm2 respectively (the 90% completeness limit, Fig. 16). The flux uncertainties are about 25-30 % and going down to fainter fluxes leads to very poor photometry.

## 9 Conclusions

Various procedures for detecting and characterizing sources were tested by simulated X-ray images. We have concentrated our attention mainly on images with XMM-Newton specific characteristics, because the problems arising from its high sensitivity and relatively large PSF are new and challenging.

We have analyzed the detection rate and the recovery of all characteristics of the input objects: flux, positional accuracy, extent measurements and the recovery of the input $`\mathrm{log}N\mathrm{log}S`$ relation. We have also investigated confusion problems in large exposures.

Concerning detection rate and characteristics reconstruction, we have shown that the VTPDETECT implementation of the Voronoi Tessellation and Percolation method is not suited to XMM-Newton images. EWAVELET provides very good detection rate and photometric reconstruction for point-like sources after a simple correction, but shows unreliable results for extended sources.

One of the best methods for point-like source detection and flux measurements is EMLDETECT but we stress again that the PSF model used for the ML procedure needs to be close to the image PSF for the most accurate photometry. Serious drawbacks are the relatively large number of spurious detections as well as the splitting of the extended sources, which we were not able to suppress even with 6 simultaneous PSF profile fits in the multi-PSF mode; this seriously hampers the analysis of the extended sources.

WAVDETECT is a flexible method giving good detections even in some complicated cases. But, here again, spurious detections are quite numerous. WAVDETECT does not assume a PSF model but requires the PSF size as a function of the encircled energy fraction and the off-axis distance in order to define the object detection cell. However, the way the detection cell is defined leads to bad photometry for extended objects.

Our choice is the MR/1+SE method. The mixed approach involving first a multiresolution iterative threshold filtering of the raw image followed by detection and analysis with SExtractor. Our tests have shown that this is the best strategy for detecting and characterizing both point-like and extended objects . Even though this mixed approach consists of two distinct steps, it is one of the fastest procedures (Tab. 9), allowing easy checks of different stages in the analysis (filtering, detection, photometry).

Without blending or confusion effects, the photometry is accurate within 10-20% for both point-like and extended objects. This uncertainty can be regarded as an intrinsic error due to the Poissonian nature of the X-ray images. For extended objects, only the MR/1+SE method gives satisfactory photometric results.

Blending between extended and point-like sources is quite serious at separations below $`30\mathrm{}`$. Better results for photometry may eventually be obtained if the intrinsic shape of the extended objects is known, and if the two objects are detected. However, in most of the cases with small separation, there is no indication of blending – which is a dangerous situation for flux reconstruction. In such cases, there may exist some spectral signatures of the effect.

The identification process of X-ray sources relies on their positional accuracy. We have shown that for objects with more than 100 counts in 10 ks exposure images and within the inner $`10\mathrm{}`$ of the field-of-view, the one sigma positional error is of the order of one half of the FWHM of the PSF ($`34\mathrm{}`$, Tab. 6). For extended objects, because of their shallower profiles and depending on the number of photons and the off-axis distance, the detected centre could even be at about $`15\mathrm{}`$ from its input position.

Comparing series of simulations with 100 ks and 10 ks in two energy bands – $`[0.52]`$ and $`[210]`$ keV, we show that the effects of confusion and completeness are absent for 10 ks, but quite significant for 100 ks in the lower energy band. Moreover, for faint fluxes, these effects tend to be masked by the large number of spurious detections with EMLDETECT. Although this method seems to give correct results for the $`\mathrm{log}N\mathrm{log}S`$ down to fainter fluxes than MR/1+SE, in real situations it is impossible to asses the contribution of the numerous spurious detections. From our simulations, we estimate that about 60-65% of the sources are lost between $`3\times 10^{16}`$ and $`6\times 10^{16}`$ erg/s/cm2 for a 100 ks exposure with the current best method (MR/1+SE).

One of the most important conclusions that will have deep cosmological impact concerns the detection and classification of extended objects. We have shown that the MR/1+SE mixed approach is capable of detecting galaxy cluster-like objects with moderate luminosity ($`L_{[210]keV}3\times 10^{44}`$ erg/s) at redshifts $`1.5<z<2`$ in 10 ks XMM-Newton simulated images. A criteria based on the half-light radius and the stellarity index classifies them correctly, with a confidence level greater than 98%.

###### Acknowledgements.

We are thankful to J.-L. Starck for many discussions regarding wavelet filtering and detections and for the MR/1 software, R. Ogley and A. Refregier for valuable comments on the manuscript, H. Bruner and J. Ballet for comments and help on XMM-SAS and EMLDETECT, E. Bertin for help on SExtractor. We thank also the referee for valuable comments and suggestions on the manuscript.


format: text

text:
# Formation of High Mass X-ray Black Hole Binaries

## 1 Introduction

Accepted wisdom in the astrophysical community (for recent relevant articles see Portegies Zwart, Verbunt, & Ergma 1997; Ergma & van den Heuvel 1998) has been that the results of burning the helium core of a star are independent of whether it is covered by a hydrogen envelope or whether it burns as a “naked” star, the envelope having been removed; e.g. by Roche Lobe overflow in main sequence or red giant stage of the star. This latter situation takes place in binary evolution where the evolution referred to in the above two articles determines the black hole progenitors to have ZAMS mass $`>40M_{}`$. Thus, this mass limit for black hole formation is adopted for the evolution of the black hole soft X-ray transients (SXTs) in these articles. On the other hand, as acknowledged in them, the observed number of SXTs requires that many, if not most, stars of ZAMS masses in the interval of $`2035M_{}`$ evolve into high-mass black holes, of mass $`\text{ }>6M_{}`$.

Brown, Lee, & Bethe (1999) offered as explanation of the large number of SXTs that stars in the mass range $`2035M_{}`$ could evolve into high-mass black holes as long as their H envelope was lifted off by the companion star only after the He core burning was complete (Case C mass transfer). If it were lifted off earlier, too small an Fe core evolved to make a high-mass black hole. In fact, Timmes, Woosley, & Weaver (1996) had found that Fe cores of Type Ib supernovae were systematically of lower mass than those of Type II, basically the same concept. The Brown, Lee, & Bethe scenario had been outlined by Brown, Weingartner, & Wijers (1996). Their results were not generally believed, however, because the calculation of Woosley, Langer, & Weaver (1995) on which they relied employed mass loss rates that were too large by a factor of 2 to 3 for the naked He stars. Thus, the present work using lower, more correct, rates is necessary.

Wellstein & Langer (2000) recalculated the evolution of naked He cores to their carbon-oxygen cores with wind loss rates lowered by a factor of 2, and in some cases by larger factors, finding that because of feedbacks, the total wind loss scaled more slowly than linearly with the multiplayer on the wind loss rate. In work that we describe here, Alexander Heger, using the Woosley & Weaver computer program Kepler, evolved the CO cores further to their Fe cores. These are not sufficiently massive to end up as high-mass black holes; in fact, the situation is nearly the same as found by Brown, Weingartner & Wijers (1996) where too large He star mass loss rates had been used. Our main objective in this article is to explain why the masses of the final Fe cores depend only insensitively on the He core mass loss rates. We interpret the Fe core mass as a good indicator of the final compact core mass because binding energy corrections and mass fallback very nearly compensate each other (Brown, Weingartner, & Wijers 1996).

The plan of this paper is as follows:

In Sec. 2 we discuss the improved measurements of wind losses as deduced from Wolf-Rayet stars.

In Sec. 3 we show that the large differences between “naked” and “clothed” He star burning can be understood in terms of the different ways in which carbon burns in the two cases.

In Sec. 4 we deal with a possible evolution of the $`10M_{}`$ black hole in Cyg X-1, suggesting that there may be an analogy with the evolution of transient sources; namely that the He core must be kept covered by H during most of its burning. This requires a very massive progenitor, a WNL star.

In Sec. 5 we give simple concluding remarks.

## 2 Wind Loss Rates in WNE’s (naked He stars)

It is the mass loss rate proposed by Langer (1989a) for hydrogen-free Wolf-Rayet stars, which replaced the previously used assumption of a constant Wolf-Rayet mass loss rate in many evolutionary calculations (Schaller et al. 1992; cf. WLW 1995), which leads to the small final masses of stars above $`35M_{}`$. This semi-empirical rate has been criticised as too high, one argument being the existence of massive black hole binaries. Indeed, recently studies of Wolf-Rayet stars show that originally measured wind losses have to be corrected downwards by a factor of 2…3 to account for their “clumpiness” (Hamann & Koesterke 1998). This is supported by polarisation measurements of the Thomson scattering, which depend linearly on the wind density (St.-Louis et al. 1993, Moffat & Robert 1994), and is in approximate agreement with the rates that would be deduced from the observed rate of increase in orbital periods for spherical mass loss

$`{\displaystyle \frac{\dot{M}}{M}}={\displaystyle \frac{\dot{P}}{2P}}`$ (1)
in Wolf-Rayet binaries. In V444 Cygni $`\dot{P}=0.202\pm 0.018\mathrm{s}\mathrm{yr}^1`$ was obtained by Khaliullin et al. (1984) and $`M_{WR}=9.3\pm 0.5M_{}`$ by Marchenko et al. (1994), resulting in

$`\dot{M}_{dyn}1\times 10^5M_{}\mathrm{yr}^1.`$ (2)
This is to be compared with the

$`\dot{M}=0.75\times 10^5M_{}\mathrm{yr}^1`$ (3)
obtained by St.-Louis et al. (1993) from the polarisation measurements. In later work Moffat & Robert (1994) arrive at a mean of $`(0.7\pm 0.1)\times 10^5M_{}\mathrm{yr}^1`$. For WNE (defined here as a Wolf-Rayet phase with vanishing surface hydrogen abundance) stars, WLW (1993) used a He wind rate of

$`\dot{M}=5\times 10^8(M_{WR}/M_{})^{2.6}\dot{M_{}}\mathrm{yr}^1`$ (4)
which would give a mass loss rate of $`1.6\times 10^5M_{}`$ yr-1 for a $`9.3M_{}`$ WR, a factor of 1.6 to 2.1 greater than eqs. (2) and (3). Given the errors in measurement, we feel safe in saying that the true wind loss rate is not less than 1/3 to 1/2 the WLW value. Furthermore, we should mention that absolute value slope of the mass loss relation Eq. (4) agrees well with the latest empirical mass loss rate determinations of Nugis & Lamers (2000).

However, decreasing the rate to 1/3 does not change our conclusions, as shown by recent calculations with such a reduced mass loss rates by Wellstein & Langer (1999), by Fryer et al. (2001), and by the calculations of Heger described later in this paper. These show that stars with initial masses below $`60M_{}`$ are unlikely to produce black holes more massive than $`10M_{}`$, from naked helium stars, although high mass black holes ($`510M_{}`$) may result from the ZAMS mass range $`2035M_{}`$ employed by Brown, Lee, & Bethe (1999) provided Case C mass transfer occurs in the binary. The relative insensitivity of the final Fe core masses to the wind-loss rate depends somewhat intricately on the way carbon burns in a naked He core, as we describe in the next section.

Before we go on to the next section we discuss, however, the possible evolution of Cyg X-1, a continuously shining black hole binary where the black hole mass is observed to be $`\text{ }>10M_{}`$. Even though massive stars are so prone to lose their hydrogen envelopes and then suffer extreme wind mass loss as hydrogen-free Wolf-Rayet stars, Langer (1987) pointed out that the most massive stars may in fact avoid this kind of evolution. The reason for doing so is that these stars, upon core hydrogen exhaustion, produce an extended region of intermediate hydrogen abundance in between the helium core and the hydrogen-rich envelope. While the latter is supposed to be lost quickly in a so called Luminous Blue Variable (LBV) stage (e.g., Stothers 2000), this intermediate region may be sufficiently massive that the star does not manage to blow it away during the ensuing so called WNL stage (defined here as a Wolf-Rayet phase with non-vanishing surface hydrogen abundance).

An example for the resulting evolutionary sequence O star $``$ LBV $``$ WNL $``$ SN is the $`100M_{}`$ sequence of Langer & El Eid (1986). Also WLW (1993) showed that the final masses of stars above $`80M_{}`$ may increase again for larger initial masses, while stars in the range $`35\mathrm{}60M_{}`$ obey the opposite trend. Even using the high mass loss rate, their $`85M_{}`$ WRB star evolved a final He core of $`9.71M_{}`$. We will estimate in Section 4 that this can be brought up to $`16M_{}`$ with the lower, more correct, mass loss rates. The reason that the He core of the WLW $`85M_{}`$ star was so much more massive than that of their $`60M_{}`$ star, was not only the higher mass of the progenitor but the fact that the WNL stage of the $`85M_{}`$ comprised $`40\%`$ of the Wolf-Rayet stage, whereas that of the $`60M_{}`$ star was $`25\%`$. During the WNL stage the star is covered by hydrogen, so that the mass loss rate is smaller than it would have been for a naked He star, and moreover the helium core mass increases with time.

Although the evolution of star as massive as $`85M_{}`$ is very uncertain, we believe that the only possibility for them to end up in $`\text{ }>10M_{}`$ black holes is an extended WNL stage, lasting a major fraction of the core helium burning. As we discuss below, these stars would have an active hydrogen burning shell source, whereas the existence of such a stage in a single massive star would probably not affect its evolution substantially, an expansion of the star due to hydrogen shell burning does give the possibility of common envelope evolution with a companion. In other words, the most massive stars are covered by H during much of their core He burning time, and this may enable them to end up as high-mass black holes.

## 3 “Naked” vs “Clothed” Helium Core Burning

### 3.1 Central Carbon Abundances vs Convective Carbon Burning

As noted earlier, Timmes et al. (1996) and Brown, Weingartner, & Wijers (1996) showed that compact cores which followed from “naked” helium stars were substantially less massive than followed from those “clothed” during helium burning with hydrogen envelopes. More recently, this has been explored in detail by Wellstein & Langer (1999) and Fryer et al. (2001). An important property of the stellar core that is significantly influenced by this difference is the carbon abundance in the core at central helium depletion. If its abundance is high enough, central carbon burning become exothermic enough to significantly overcome the neutrino losses and burn in a convective core. Otherwise the first convective burning of carbon can only appear in a shell source. The carbon abundance also has significant effect on the location, duration, and extent of the carbon burning shells. For example, as the carbon abundance goes to zero, no convective shell burning occurs either. The lower the carbon abundance, the further out the first shell forms. The location of the last carbon-burning shell sets the size of the carbon-free core which determines the final evolution of the star and the sizes of the iron core (possible direct black hole formation) and the silicon core (possible black hole formation by fall-back).

There are two major factors that contribute to the resulting carbon abundance after core helium depletion. The first is the dependence on the mass of the helium core and results chiefly from the different behavior of the 12C$`(\alpha ,\gamma )^{16}`$O reaction, as described clearly by Weaver & Woosley (1993). Since this is important to our development, we repeat the argument, supplemented by more recent developments.

Carbon is formed in the triple $`\alpha `$-process; since this is a three-body process, it depends on the density as $`\rho ^2`$. Carbon is removed, when possible, by the 12C$`(\alpha ,\gamma )^{16}`$O process, which goes linearly with $`\rho `$. Now, for massive helium cores the central density scales roughly as $`M^{1/6}`$. With increasing helium core mass this implies that carbon formation, which goes as $`\rho ^2`$, will be cut down, compared with 12C$`(\alpha ,\gamma )^{16}`$O, which goes as $`\rho `$. Here, we also note that both the 12C$`(\alpha ,\gamma )^{16}`$O reaction and the triple $`\alpha `$-process have very similar temperature dependence in the temperature regime of central helium burning, $`2\times 10^8`$K, i.e., the carbon production does not depend on temperature in this regime.

This means that there is a mass above which the post-helium burning central carbon abundance is low enough to skip central carbon burning. Woosley & Weaver (1995) found this transition at a mass of 19$`M_{}`$. However, this limit strongly depends on the still uncertain 12C$`(\alpha ,\gamma )^{16}`$O reaction rate and the stellar evolution model, in particular the prescription used for convection or mixing processes in general, and also mass loss. On this last issue we elaborate later in more detail.

A bare helium star behaves differently from a clothed one in that its helium core does not grow in mass due to hydrogen shell burning, as is the case in clothed stars, but rather shrinks due to mass loss from the surface – the whole star is only a helium core, a Wolf-Rayet star. Generally, the convective core tends to comprise a larger fraction of the helium core, by mass, as the central helium abundance decreases, due to the density-dependence outlined above. This statement depends, again, somewhat on the description of semiconvection, but probably applies to most models with either Schwarzschild convection (e.g., Schaller et al. 1992), overshooting or fast semiconvection (e.g., Woosley & Weaver 1995), or rotationally induced mixing (Heger, Langer, Woosley 2000). Additionally, the mass fraction of the convective core also increases with the mass of helium core. However, bare helium cores (early-type Wolf-Rayet stars) experience mass loss rates that are sufficiently strong that the convective core actually tends to shrink in the run of its evolution, in particular also towards the end of helium burning, rather than grow, as it does in a clothed star. This is despite the fact that it still comprises an increasingly larger fraction of the remaining total helium core. The important point is that a growth in mass of the convective core injects new helium into this convective core, while when the mass is constant or decreases, this injection does not occur.

As the triple alpha reaction depends on the third power of the helium mass fraction it loses against the 12C$`(\alpha ,\gamma )^{16}`$O reaction toward the end of central helium burning, i.e., carbon is mostly burned rather then produced toward the end of central helium burning. That switch typically appears at a central helium mass fraction of around 10–20 %. Most importantly, as can be seen from the central carbon abundances at the end of He burning in Table 1, the He fraction is too low to burn the final carbon, which remains at a central abundance $`\text{ }>0.30`$ for all mass stars. This is roughly double the carbon abundance necessary for convective carbon burning. It is seen in Table 2 that the central carbon abundance in naked He stars goes down only slowly with decreasing mass loss rate, although by the time the mass loss is reduced by a factor of 6 the convective burning stops. In the clothed stars the growth of the core and its accompanied injection of helium after this time thus leads to a further decrease of carbon as compared to the bare helium cores that do not have this additional supply of helium.

In Fig. 1 we show calculations carried out by Tom Weaver (priv. com. 1995) using the KEPLER code (Weaver, Zimmerman, & Woosley 1978; Woosley & Weaver 1995) on the central carbon abundance after core helium depletion for “clothed” (single) stars without (any) mass loss as function of ZAMS mass. The rapid drop in $`C_c`$ at ZAMS mass $`M_{ZAMS}20M_{}`$ causes the disappearance of convective carbon burning.

### 3.2 Iron Core Masses

It is generally known that the central entropy rises with the mass of the star; so that the more entropy is carried off, the less massive will be the final Fe core. This statement can be made more quantitatively by using the Bethe, Brown, Applegate, & Lattimer (1979) conclusion that the entropy per nucleon in the Fe core is $`\text{ }<1`$, in units of the Boltzmann constant. For a ZAMS $`20M_{}`$, the central entropy in our Fe core is 0.76 per nucleon when the He core is burned clothed, 0.775 when burned naked; i.e., essentially the same. Thus, the total entropy of the Fe core is $`N_N`$ considering some gradient towards higher entropy at the outer layers of the core, where $`N_N`$ is the number of nucleons in it. In this way one can see directly that the entropy that is carried off during the evolution of the Fe core will diminish its mass.

In Fig. 2 we show iron core masses at the time of iron core collapse for a finely spaced grid of stellar masses (Heger, Woosley, Martinez-Pinedo, & Langanke 2001). Filled circles and crosses correspond to the core masses of “clothed” single stars. The circles were calculated by Tom Weaver (1995, priv. com.) using the same physics as in Woosley & Weaver (1995), whereas the crosses employ the improved weak rates by Langanke & Martinez-Pínedo (2000) for electron capture and beta decay. Although the latter are much smaller than the rates by Fuller, Fowler, & Newman (1985) used by Woosley & Weaver (1995), the final Fe core masses are not much changed, for reasons described by Heger, Woosley, Martinez-Pinedo, & Langanke (2001).

For “clothed” single stars, one sees an anticorrelation between the peak in Fe core masses at ZAMS mass $`23M_{}`$ and the minimum in $`C_c`$ at ZAMS mass $`21M_{}`$. In other words, the peak around $`23M_{}`$ in ZAMS masses occurs just where the central carbon abundance at the end of helium core burning is at its minimum of $`10\%`$, too low for convective core burning or the formation of a carbon-burning convective shell close to the center of the star.

When sufficient carbon remains, it will be burned convectively in reactions such as 12C$`+^{12}`$C$`^{24}`$Mg, 20Ne$`+\alpha `$ etc., at a temperature of $`7080`$ keV, several times higher than the $`20`$ keV needed for 12C$`(\alpha ,\gamma )^{16}`$O. If the carbon abundance is lower, the carbon burning phases are reduced and typically bigger carbon-free cores result for otherwise similar stars. Therefore less energy is carried away by neutrinos and the entropy in the core stays higher. More massive silicon and iron cores can form in this case (Boyes, Heger, and Woosley 2001). This is clearly seen in Figures 1 & 2.

Although the central carbon abundance goes up again for ZAMS mass $`25M_{}`$, and the Fe core mass goes down, by such high masses the gravitational energy of the stellar envelope is large and it is difficult for the shock energy after collapse to blow it off, so we believe that these will go into high-mass black holes, as well as the stars of ZAMS mass $`2023M_{}`$.

### 3.3 Formation of High-Mass Black Holes

Although there is a lot of uncertainty in the literature about the ZAMS masses which end up in high mass black holes, Bethe & Brown (1999) argued that this should occur at the proto compact core mass of

$`M_{PC}1.8M_{}.`$ (5)
Further uncertainty comes in relating the Fe core mass to the compact core mass. It can, however, be seen from Table 3 of Brown, Weingartner, & Wijers (1996) that the estimated fallback mass following separation in the supernova explosion roughly compensates for the binding energy increase in going from the Fe core to the compact object. Thus, we will assume the compact object masses to be the same as those of the Fe cores. In addition to the uncertainty in this relation, there is the further uncertainty that the fallback mass, which is dependent on the density structure in the layers above the Fe core (Fryer 1999; Janka 2001), is unknown, since consistent supernova explosions have not been achieved.

None the less, we believe it useful to move ahead with our estimated $`M_{PC}`$ which indicates that single stars with ZAMS masses $`\text{ }>20M_{}`$ go into high-mass black holes. Brown, Lee, & Tauris (2001) recently showed that given presently accepted wind losses (Schaller et al. 1992) the high mass black holes in the transient sources with main sequence companion can only be evolved in the region of ZAMS mass $`20M_{}`$, so our choice of $`M_{PC}`$ is supported to this extent by evolutionary arguments.

We note briefly that SN 1987A, if it formed a black hole, went into a low mass black hole of mass $`1.5M_{}<M_{BH}<1.8M_{}`$ according to the Brown & Bethe (1994) estimates. The black hole would have a low mass in this interval because the He envelope was blown off in the delayed explosion. Note from Fig. 2 of main sequence masses that there is only a narrow interval from $`18M_{}`$ to $`20M_{}`$ in which this could happen given the above interval.

Schaller et al. (1992) use a 12C$`(\alpha ,\gamma )^{16}`$O S-factor of S(300 keV) $`100`$ keV barns as compared with the Woosley & Weaver 170 keV barns. They bring their central abundance down to 0.16 at the end of core He burning only for a ZAMS $`25M_{}`$ star, so presumably this would be the mass at which the Fe cores would begin to rise rapidly in mass with further evolution. We believe the Weaver & Woosley value for the 12C$`(\alpha ,\gamma )^{16}`$O rate to be more correct (Boyes, Heger, & Woosley 2001), placing the narrow interval of ZAMS masses from which SN 1987A can be evolved correctly. In fact, most recently experiments including both E1 and E2 components have been carried out (Kunz et al. 2001), obtaining $`S_{tot}^{300}=(165\pm 50)`$ keV barns.

## 4 Possible Evolution of Cyg X-1

### 4.1 Space Velocity vs the Mass Loss in the Formation of Cyg X-1

According to the calculations of Wellstein & Langer (1999) and Fryer et al. (2001), it is unlikely that a black hole as massive as the $`10M_{}`$ core (Herrero et al. 1995) in Cyg X-1 can be evolved from a naked He star (see Table 1 and Fig. 2). Although the calculations of naked He stars were extended up to only ZAMS $`60M_{}`$ stars, the He winds scale with the 2.6 power of the mass, so that higher mass He stars would be expected to behave in the same way. In fact, Wellstein & Langer (1999) and Fryer et al. (2001) were unable to evolve a $`10M_{}`$ black hole, whereas the high space velocity of Cyg X-1, $`49\pm 14`$ km s-1 (Kaper et al. 1999) indicates substantial mass loss in a Blaauw-Boersma kick (Blaauw 1961; Boersma 1961) in the black hole formation, increasing the necessary He core mass. We now estimate this mass loss, following the development of Nelemans et al. (1999). It is reasonable to evolve Cyg X-1 analogously to the soft X-ray black hole transient sources, the difference being in the high mass companion star which presumably makes Cyg X-1 shine continuously.

The high space velocity of Cyg X-1 can be explained by mass loss in the black hole formation, because this loss from the black hole is somewhat off from the center of gravity of the system. From Nelemans et al. (1999), the runaway velocities from symmetric SNe (See Appendix for the derivation)

$`v_{sys}`$ $`=`$ $`213\times \mathrm{\Delta }M\times m\times P_{recirc}^{1/3}\times (M_{BH}+m)^{5/3}\mathrm{km}\mathrm{s}^1`$ (6)
where masses are in $`M_{}`$, $`P`$ in days, and $`\mathrm{\Delta }M=M_{He}M_{BH}`$ with $`M_{He}`$ the He star mass of the Black hole progenitor. Here $`P_{recirc}`$ is the re-circularized orbital period after Blaauw-Boersma kick, and we assume no orbital evolution ($`P_{recirc}=P_{obs}`$) since the beginning of the mass transfer phase, and neglect small eccentricity as in Nelemans et al. (1999). By putting in average values,

$`v_{sys}`$ $`=`$ $`8.32\times \mathrm{\Delta }M\mathrm{km}\mathrm{s}^1.`$ (7)
In order to obtain the observed velocity, we need $`\mathrm{\Delta }M5.9M_{}`$, indicating the black hole progenitor mass to be $`M_{He}16M_{}`$. Note that this $`\mathrm{\Delta }M`$ and $`M_{He}`$ are in the same ballpark as those estimated by Nelemans et al. for Nova Scorpii 1994(GRO J1655$``$40). This suggests that the evolution of Cyg X-1 is similar to that of Nova Scorpii 1994, except that the progenitor mass of the black hole must be much higher in Cyg X-1 and the companion is an O-star rather than an F-star.

The current orbital separation of Cyg X-1 is

$`a_{\mathrm{now}}4.2R_{}\left({\displaystyle \frac{P_{orb}}{\mathrm{days}}}\right)^{2/3}\left({\displaystyle \frac{M_{BH}+M_O}{M_{}}}\right)^{1/3}40R_{}.`$ (8)
Before the explosion, the orbital separation of the black hole progenitor and companion star was

$`a_{\mathrm{preSN}}{\displaystyle \frac{a_{\mathrm{now}}}{1+e_2}}33R_{}`$ (9)
where $`e_2=\mathrm{\Delta }M/(M_{\mathrm{BH}}+M_O)0.21`$ is the eccentricity right after the supernova explosion. For the derivation, see Appendix. This is only slightly larger than the $`30R_{}`$ estimated by Bethe & Brown (1999) for the minimum initial separation of the O-star progenitors. If the progenitors were initially closer together they would merge already at this stage of evolution.

As noted by Bethe & Brown (1999), the two O-star progenitors must have initially had a separation of nearly $`33R_{}`$. Substantial mass loss, especially by the WR star would have been expected to widen the orbit substantially

$`{\displaystyle \frac{a_f}{a_i}}={\displaystyle \frac{M_i}{M_f}}23`$ (10)
where $`M_i`$ and $`M_f`$ are the initial and final system masses. In order to tighten the orbit to the $`33R_{}`$ we find before the massive star goes into a black hole, there must have been some period of nonconservative mass transfer.

### 4.2 Extended WNL Stage of Black Hole Progenitor

Langer (1989b) finds that towards the end of core hydrogen burning the radius of a massive star may be much larger than its radius on the zero age main sequence. For example, a $`100M_{}`$ star increases its radius from $`13R_{}`$ at H ignition by a factor $`4`$ to $`53R_{}`$ at core H exhaustion. In the evolutionary phase between central H- and He-burning massive stars develop an intermediate fully convective zone just above the H-burning shell which reaches its maximal spatial extent and $`18M_{}`$ in mass for the $`100M_{}`$ star (Langer 1987). There may be further expansion in a hydrogen shell burning stage unless large mass loss sets in. But the latter may begin only later in the LBV stage. In any case, there will be a period before or in the LBV stage where the hydrogen is being propelled outwards, but does not yet have enough velocity to escape. A companion to the very massive star will, if in this region, couple hydrodynamically to the hydrogen in the manner explicitly worked out by Bethe & Brown (1999), and furnish energy to it from a drop in its gravitational binding; i.e., through some mechanism resembling common envelope evolution. We cannot be more explicit, because of the uncertainty in the post main sequence evolution of the very massive stars, but some nonconservative mass transfer seems to be necessary in the evolution of a binary such as Cyg X-1. We admit to being in somewhat of a predicament, however, since we do not want a common envelope evolution which removes the hydrogen before the helium burning is well underway.

It is clear that the black hole in Cyg X-1 cannot have evolved from a naked He star. An important clue to its possible evolution is offered by the evolution of WLW (1993) of a ZAMS $`85M_{}`$ star, where the WR stage began with a WNL stage which took $`40\%`$ of the total WR time. (Their $`60M_{}`$ star also had a WNL stage, but it took only $`25\%`$ of the WR time.) During this time the WR is covered by hydrogen and the attendant mass loss rate is much lower than that for a naked He star.

Using the times for each WR stage as WLW 1995, we calculate in Table 3 the masses for ZAMS $`60M_{}`$ and $`85M_{}`$ stars. Of course there is feedback on these masses from the change in times of each stage with altered mass loss rate. In order to take the feedback into account, we use the ratios of $`M_{CO}`$ cores calculated by Fryer et al. (2001) for a ZAMS $`60M_{}`$ star for various mass loss rates to the original WLW (1995) one as plotted in Fig 3.

We note from Table 3 that a ZAMS $`60M_{}`$ star with WR stage mass loss rate reduced by 1/3 loses $`1.7M_{}`$ during its WNL stage, at an average rate of $`1.5\times 10^5M_{}`$ yr-1. This is an order of magnitude less than it would lose as a WNE (bare) He star, and explain why the final He star mass can be as large as $`13.5M_{}`$. The WNL mass loss rate is $`2.25`$ times the main sequence mass rate for a ZAMS $`60M_{}`$ star of Castor, Abbott, & Klein (1975).

In Fig. 3 we plot our ratios from Table 3 for He star masses alongside the Fryer et al. (2001) ratios of CO core masses. The divergence in our favored range of $`1/21/3`$ is not large, so we believe our procedure in calculating masses is justified. We see from Table 3 that we need a mass loss rate of $`0.4`$ times the WLW (1995) one in order to lose $`6M_{}`$ in the explosion and to be left with an $`10M_{}`$ black hole.

### 4.3 Discussion on Cyg X-1 Type Objects

We want to point out that there is a strong metallicity dependance in our model. As the winds of WNL stars are likely radiation driven (Hamann et al. 2000), we can expect them to weaken with smaller metallicity. Thus, the lower limit of the initial mass range from which on the the final stellar mass can grow again will be decreasing. Therefore, we expect more Cyg X-1 type systems at lower metallicity. We also would like to point out that there is dependance in mixing processes, in particular rotation, which can increase the lower limit of the initial mass range.

Using standard population synthesis, Bethe & Brown (1999) estimated that there should be $`7`$ Cyg X-1 like objects in the Galaxy, where we see only one (leaving out LMC X-1 and LMC X-3). Although the $`\gamma `$-rays from such an object easily penetrate the Galactic disc, they would not be distinguished from those from gamma ray bursts (Maarten Schmidt, private communication 2001).

As possible WR progenitor of the Cyg X-1 black hole, we suggest “WR22”; the most massive Wolf-Rayet ever weighed (Rauw et al. 1996). The minimum Wolf-Rayet mass is $`72M_{}`$ and the mass ratio to O-star is 2.78, giving an O-star mass of $`>26M_{}`$. The spectrum of the WN7 exhibits a considerable amount of hydrogen ($`X_H40\%`$, Hamann et al. 1991, Crowther et al. 1995), strong WR emission lines and absorption lines that belong to the WR component (Niemela 1973; Moffat & Seggewiss 1978). Rauw et al. (1996) identify the companion as an O-star. With $`80`$ day period, the binary is, however, much too wide in order to narrow sufficiently in common envelope evolution to produce a Cyg X-1 type binary, according to the Bethe & Brown (1999) estimates.

## 5 Conclusion

Our conclusion to all of the above is simple: In order to evolve high-mass black holes at solar metallicity, the He core of the massive star must be covered by hydrogen at least most of the time, while it burns.

For SXTs, the stars in the mass range $`2035M_{}`$ could evolve into high-mass black holes with the condition that their H envelope was lifted off by the companion star only after the He core burning was complete. For more massive stars like Cyg X-1 type objects, an extended WNL stage could drive the formation of high-mass black holes.

## Acknowledgments

We thank Stan Woosley for many helpful discussions and are indebted to Tom Weaver for supplying us with his grid of stellar evolution models. We wish to thank Maarten Schmidt for communication about the indistinguishability of further Cyg X-1 like objects from gamma ray bursts. This work was partially supported by the U.S. Department of Energy under Grant No. DE–FG02–88ER40388, by the Alexander von Humboldt-Stiftung through Grant FLF-1065004, by 2000-2001 KIAS Research Fund, and by the Deutsche Forschungsgemeinschaft through Grants La 587/15 and La 587/16.

## Appendix. Blaauw-Boersma Kick

In this Appendix, the details of the Blaauw-Boersma kick (Blaauw 1961, Boersma 1961) are summarized. Before the explosion,111Subscript 1 (2) indicates the properties before (right after) the explosion, and subscript “now” indicates the current properties observed. from Fig. 4 (a), $`R_{He}`$, $`v_{He}`$, and $`\omega _1`$ are given as

$`R_{He}`$ $`=`$ $`\left({\displaystyle \frac{m}{M_{He}+m}}\right)a_1`$
$`v_{He}`$ $`=`$ $`{\displaystyle \frac{1}{a_1}}\left(GmR_{He}\right)^{1/2}`$
$`\omega _1`$ $`=`$ $`\left({\displaystyle \frac{m}{M_{He}+m}}\right)\left({\displaystyle \frac{Gm}{R_{He}^3}}\right)^{1/2}.`$ (11)
After the explosion, we assume that $`\mathrm{\Delta }MM_{He}M_{BH}`$ was lost from the BH progenitor without interacting with the binary system. The momentum lost from the binary system should be compensated by the space velocity of the new binary system.

$`\mathrm{\Delta }p=\mathrm{\Delta }Mv_{He}=(M_{BH}+m)v_{sys}`$ (12)
which gives the space velocity of the new binary system (c.m. motion of B)

$`v_{sys}=\left({\displaystyle \frac{\mathrm{\Delta }M}{M_{BH}+m}}\right)v_{He}.`$ (13)
Right after the explosion, in the new c.m. frame of BH and O-star binary, by defining the semi-major axis $`d_2`$ and the eccentricity $`e_2`$, we have

$`\omega _2`$ $`=`$ $`\left({\displaystyle \frac{m}{M_{BH}+m}}\right)\left({\displaystyle \frac{Gm}{d_2^3}}\right)^{1/2}`$
$`R_{BH}`$ $`=`$ $`\left({\displaystyle \frac{m}{M_{BH}+m}}\right)a_1`$
$`v_{BH}`$ $`=`$ $`v_{He}+v_{sys}=\left({\displaystyle \frac{M_{He}+m}{M_{BH}+m}}\right)v_{He}.`$ (14)
Since the $`v_{BH}`$ is greater than that of the circular motion with radius $`R_{BH}`$

$`v_{BH}`$ $`=`$ $`\left({\displaystyle \frac{M_{He}+m}{M_{BH}+m}}\right)^{1/2}v_{circ}>v_{circ}`$
$`v_{circ}`$ $`=`$ $`\left({\displaystyle \frac{Gm^2}{a_1(M_{BH}+m)}}\right)^{1/2},`$ (15)
one can conclude that $`R_{BH}`$ is at the minimum of the orbital separation, i.e.

$`R_{BH}=d_2(1e_2).`$ (16)
Note that $`\stackrel{}{v}_{BH}`$ and $`\stackrel{}{R}_{BH}`$ are perpendicular right after the explosion as in Fig. 4 (c).

The orbital angular momentum of the BH can be expressed using $`a_2`$ and $`e_2`$

$`l_a`$ $`=`$ $`\left[Gm\mu ^2d_2(1e_2^2)\right]^{1/2}=\left[Gm\mu ^2R_{BH}(1+e_2)\right]^{1/2}`$ (17)
or

$`l_b`$ $`=`$ $`M_{BH}R_{BH}v_{BH}=\mu a_1v_{BH}=\mu \left(GmR_{BH}\right)^{1/2}\left({\displaystyle \frac{M_{He}+m}{M_{BH}+m}}\right)^{1/2}`$ (18)
where $`\mu `$ is the reduced mass

$`\mu ={\displaystyle \frac{M_{BH}m}{M_{BH}+m}}.`$ (19)
By equating $`l_a=l_b`$ we have

$`e_2={\displaystyle \frac{\mathrm{\Delta }M}{M_{BH}+m}}.`$ (20)
We assume that there is no angular momentum lost during the circularization process by tidal locking after SN explosion. From angular momentum conservation

$`l_{now}`$ $`=`$ $`\left(Gm\mu ^2d_{now}\right)^{1/2}=\left[Gm\mu ^2d_2(1e_2^2)\right]^{1/2}.`$ (21)
Therefore, from

$`d_{now}`$ $`=`$ $`\left({\displaystyle \frac{m}{M_{BH}+m}}\right)a_{now}`$ (22)
$`=`$ $`d_2(1e_2^2)=R_{BH}(1+e_2)=\left({\displaystyle \frac{m}{M_{BH}+m}}\right)a_1(1+e_2),`$
we have the relation of the orbital separations

$`a_{now}=a_1(1+e_2).`$ (23)
Using this relation, one have

$`\omega _1`$ $`=`$ $`(1+e_2)^2\omega _{now}`$
$`\omega _2`$ $`=`$ $`(1e_2^2)^{3/2}\omega _{now}`$ (24)
or the period relations

$`P_{now}=(1e_2^2)^{3/2}P_2=(1+e_2)^2P_1.`$ (25)
Now Eq. (13) can be expressed in terms of observed quanties using

$`a_{now}`$ $`=`$ $`\left({\displaystyle \frac{P_{now}}{2\pi }}\right)^{2/3}\left[G(M_{BH}+m)\right]^{1/3}`$
$`v_{He}`$ $`=`$ $`\left({\displaystyle \frac{1+e_2}{a_{now}}}\right)^{1/2}\left({\displaystyle \frac{Gm^2}{M_{He}+m}}\right)^{1/2}=\left[{\displaystyle \frac{1}{a_{now}}}\left({\displaystyle \frac{Gm^2}{M_{BH}+m}}\right)\right]^{1/2}`$ (26)
$`=`$ $`{\displaystyle \frac{(2\pi G)^{1/3}m}{(M_{BH}+m)^{2/3}P_{now}^{1/3}}}.`$
Finally the Blaauw-Boersma kick velocity becomes

$`v_{sys}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }M}{M_{BH}+m}}v_{He}={\displaystyle \frac{(2\pi G)^{1/3}\mathrm{\Delta }Mm}{(M_{BH}+m)^{5/3}P_{now}^{1/3}}}`$ (27)
$`=`$ $`213\left({\displaystyle \frac{\mathrm{\Delta }M}{M_{}}}\right)\left({\displaystyle \frac{m}{M_{}}}\right)\left({\displaystyle \frac{P_{now}}{\mathrm{days}}}\right)^{1/3}\left({\displaystyle \frac{M_{BH}+m}{M_{}}}\right)^{5/3}\mathrm{km}\mathrm{s}^1.`$