# Distribution of compact object mergers around galaxies

## 1 Introduction

The recent discoveries of X-ray afterglows of gamma-ray bursts by the Beppo SAX satellite \[Costa et al., 1997\] have revolutionized the approach to these phenomena. For the first time since their discovery 30 years ago \[Klebesadel et al., 1968\] gamma ray bursts have been identified with sources at other wavelengths. In consequence optical afterglows have been discovered \[Groot et al., 1997a\], which lead to identification of host galaxies \[Groot et al., 1997b\]. At the time of writing more than a dozen afterglows have been identified. The optical lightcurves of the GRB afteglows decay as a power law $`t^\alpha `$ with the typical values of the index $`\alpha =1.1`$ to $`1.3`$. In some cases the host galaxies have been found by observing the flattening of of the lightcurve. The underlying steady flux is identified as the emission of the host galaxy. In a few cases the host galaxies themselves were found as extended objects. This allowed to measure the offset between the location of the gamma-ray burst and the center of the host galaxy. A list five GRBs and the offsets from their host galaxies is shown in Table 1. The offsets are generally small and the afterglows lie directly on the host galaxies. In other cases where the host galaxy have been only found from the flattening of the lightcurve we know that the location of the host galaxy and the GRB do not differ substantially from the simple fact that the host galaxies have been identified.

There are two basic categories of theoretical models of the central engines of gamma-ray bursts within the cosmological model. The first class connects gamma-ray bursts with mergers of compact objects, e.g. neutron stars and/or black holes. There exist numerous scenarios in this set of models, some of them link GRBs with the coalescence of a black hole neutron star binary \[Narayan et al., 1992\]. In other models like in the recent paper by \[Kluźniak and Ruderman, 1998\] the GRB events are related to mergers of two neutron stars. Compact object merger model provides enough energy to power a GRB, and it has been showed in the numerical simulations \[Kluźniak and Lee, 1998\] that the coalescence may last up to a second. The analytical estimations of \[Portegies Zwart, 1998\] extend this timescale up to a minute. The second class of models \[Paczynski, 1998\] relates GRBs to explosions of supermassive stars, so called hypernovae. A direct prediction in this class of models is that gamma-ray bursts are related to the star forming regions.

The relation between the GRB location and the host galaxies does not have to be true in the models that relate GRBs to compact object mergers. Tutukov and Yungelson \[Tutukov and Yungelson, 1993\] have calculated the compact object merger rates, and also found that compact object binaries may travel the distances up to $`1000`$kpc before merging. In a more detailed study \[Bloom et al., 1998a\] calculated a population of compact object binaries using the population synthesis method of Pols and Marinus \[Pols and Marinus, 1994\] and then calculated the spatial distribution of mergers in the potentials of galaxies for w few representative masses. They found that approximately $`15`$% of mergers take place outside the host galaxies). They have used the a Maxwellian kick velocity distribution with $`\sigma _v=190`$ km s<sup>-1</sup> \[Hansen and Phinney, 1997\]

In this work we use the the population synthesis code based on \[Bethe and Brown, 1998\] and extended by Belczyński and Bulik \[Belczynski and Bulik, 1999\]. We concentrate on the dependence of the properties of the compact object binaries on the parameters used in the population synthesis code. We find that the most important parameter that determines the population of compact object binaries is the kick velocity a neutron star receives at birth, however this distribution is poorly known. Iben and Tutukov \[Iben and Tutukov, 1996\] claim that the properties of pulsars can be explained by only the recoil velocities with no need for the kicks. Blaauw and Ramachandran \[Blaauw and Ramachandran, 1998\] find that a single kick velocity of $`200`$km s<sup>-1</sup> suffices to reproduce the pulsar population. Cordes and Chernoff \[Cordes and Chernoff, 1997\] proposed a weighted sum two Gaussians: 80 percent with the width $`175`$ km s<sup>-1</sup> and 20 percent with the width $`700`$ km s<sup>-1</sup>.

We outline the model of the binary evolution and propagation in a galactic potential in section 2. The results of the calculation are presented in section 3 and we discuss them in section 4.

## 2 Model

### 2.1 Binary evolution

In order to study the spatial distribution of compact object mergers we use the population synthesis method. We use the population synthesis code \[Belczynski and Bulik, 1999\] which concentrates on the population of massive star binaries, i.e. those that may eventually lead to formation of compact objects and compact object binaries. We include the evolution of the binaries due to interaction and mass transfer and also the kicks that a newly born neutron star receives in supernova explosion. A binary may be disrupted in each of the supernova events. The surviving binaries obtain center of mass velocities, which change their trajectories and may even eject them from their galaxy.

While the evolution of single stars depends only on their mass and metallicity the evolution of binaries is also a function the initial orbit (semimajor axis $`a`$, and eccentricity $`e`$) of the two stars. We assume that the distribution of the initial parameters can be expressed as a product of distributions of four parameters: the larger star (primary) mass $`M`$, the mass ratio of the less massive to the more massive star in the binary $`q`$, and the orbital parameters $`a`$ and $`e`$, i.e that this quantities are independent. The distribution of primary masses used here is \[Bethe and Brown, 1998\]

$$\mathrm{\Psi }(M)M^{3/2},$$

and we adopt a flat distribution of the mass ratio $`q`$. The semi major axis distribution is scale invariant, i.e.

$$\mathrm{\Gamma }(a)a^1$$

with the limits $`6R_{}<a<6000R_{}`$, and we draw the eccentricity from a distribution $`\mathrm{\Xi }(e)=2e`$.

We assume that the kick velocity distribution is a three dimensional Gaussian, and parameterize it with its width $`\sigma _v`$, i.e.

$$p(v)=\frac{4}{\sqrt{\pi }}\sigma _v^3v^2\mathrm{exp}\left(\frac{v^2}{\sigma _v^2}\right).$$
(1)   
We generate population of compact object binaries for a few values of $`\sigma _v`$ in order to asses the sensitivity of our results to this parameter.

We describe the mass transfer in the common envelope evolution by the common envelope parameter $`\alpha _{CE}`$ (see e.g. \[Vrancken et al., 1991\]), and we use an intermediate value of $`0.8`$ for this parameter. In this type of evolution the more massive star looses its envelope and becomes a helium star with mass approximately 30% of its initial value. The $`\beta `$ parameter which describes the specific angular momentum of the material expelled from the binary in the Roche lobe overflow phase is set to $`\beta =6`$ \[Pols and Marinus, 1994\]. Accretion onto a neutron star in a binary is treated as Bondi-Hoyle accretion and we use the formalism developed by \[Bethe and Brown, 1998\] to find the amount of mass accreted onto the neutron star, and the final orbital separation. Systems with nearly equal masses evolve at the similar speed, and loose the common envelope, shrinking their orbit at the same time. For a more detailed description of the population synthesis code see \[Belczynski and Bulik, 1999\].

We assume that a neutron star with mass of $`1.4M_{}`$ is formed in each supernova explosion. We draw a random time in the orbital motion to obtain the position on the orbit when the supernova explodes. The remaining mass of the envelope is ejected from the system, and the newly formed neutron star receives a kick.We verify whether the system is still bound after the explosion. For bound systems we find the parameters of the new orbit and the kick velocity the whole binary receives

$$\mathrm{\Delta }\stackrel{}{V}=\frac{M_2^iM_2^f}{M_1+M_2^f}(\stackrel{}{v}_2^i+\stackrel{}{v}_{kick})$$

where $`M_1`$ is the mass of the companion, $`M_2^i,M_2^f`$ are the initial and final masses of the supernova, $`\stackrel{}{v}_2`$ is the orbital velocity of the supernova at the time of explosion. After each supernova expolsion we verify whther the system survives as a binary.

A compact object binary loses its energy through gravitational radiation. The time to merge is \[Peters, 1964\]

$$t_{mrg}=\frac{5c^5a^4(1e^2)^{7/2}}{256G^3Mm(M+m)}\left(1+\frac{73}{24}e^2+\frac{37}{96}e^4\right)^1,$$
(2)   
where $`a`$ is the semi major axis of the orbit, $`e`$ is its eccentricity, and $`M,m`$ are the masses of the compact objects.

### 2.2 Orbit in Potentials of Galaxies

Since little is known about the host galaxies of gamma-ray bursts, in particular of their types and masses, we will present two extreme cases: (i) propagation in a potential of large spiral galaxy like the Milky Way, and (ii) propagation in empty space, corresponding to GRBs originating e.g. in globular clusters. In the latter case we assume that all binaries originate in one point, and travel due to kicks described above and there is no gravitational potential.

The potential of a spiral galaxy can be described as a sum of three components: bulge, disk, and dark matter halo. A convenient way to describe the Galactic potential has been proposed by \[Miyamoto and Nagai, 1975\], while a series of more detailed models were constructed by \[Kuijken and Gilmore, 1989\] and used in modeling the galactic halo population of neutron stars \[Bulik and Lamb, 1995, Bulik et al., 1998\]. The \[Miyamoto and Nagai, 1975\] potential for a galactic disk and bulge is

$$\mathrm{\Phi }(R,z)=\frac{GM}{\sqrt{R^2+(a_i+\sqrt{z^2+b_i^2})^2}}$$

where $`a_i`$ and $`b_i`$ are the parameters, $`M`$ is the mass, and $`R=\sqrt{x^2=y^2}`$. The dark matter halo potential is spherically symmetric

$$\mathrm{\Phi }(r)=\frac{GM_h}{r_c}\left[\frac{1}{2}\mathrm{ln}\left(1+\frac{r^2}{r_c^2}\right)+\frac{r_c}{r}\mathrm{atan}\left(\frac{r}{r_c}\right)\right]$$

corresponds to a mass distribution $`\rho =\rho _c/[1+(r/r_c)^2]`$. The mass of such halo is infinite, so we introduce a cutoff value $`r_{cut}=100`$kpc above which the density of the halo falls to zero. While the details of the model of galactic potential are not important for this study we have to adopt a particular value of the masses and sizes of each of the components. We use the values of the parameters as determined by \[Blaes and Rajagopal, 1991\] for the Milky Way: $`a_1=0`$kpc, $`b_1=0.277`$kpc, $`a_2=4.2`$kpc, $`b_2=0.198`$kpc, $`M_1=1.12\times 10^{10}M_{}`$, $`M_2=8.78\times 10^{10}M_{}`$, $`r_c=6.0`$kpc, and $`M_h=5.0\times 10^{10}M_{}`$.

We assume that the distribution of binaries in our model galaxy follows the mass distribution in the young disk \[Paczyński, 1990\], that is

$$P(R,z)dRdz=P(R)dRp(z)d(z).$$

The radial distribution is exponential with

$$P(R)e^{R/R_{exp}}R$$

with $`R_{exp}=4.5`$kpc and extends to $`R_{max}=20`$kpc. The vertical distribution is $`p(z)e^{z/z_{exp}}`$ and $`z_{exp}=75`$pc. We note that this is not a self consistent approach: the density inferred from the disk potential is not the same as the density of binaries. However in this work we are not interested in determining high accuracy positions around the host galaxy, and rather with an estimate of the general properties of the distribution of compact object mergers.

Each binary moves initially with the local rotational velocity in the galactic disk. After each supernova explosion we add an appropriate velocity, provided that the system survives the explosion. We calculate the orbit of each system until it merger time provided that the merger time is smaller than the Hubble time (15 Gyrs here).

## 3 Results

The kick velocity distribution is not very well known. Therefore, we use the population synthesis code with four values of the kick velocity distribution width: with no kick velocities $`\sigma _v=0`$km s<sup>-1</sup>, and with $`\sigma _v=200,\mathrm{\hspace{0.17em}400},\mathrm{\hspace{0.17em}800}`$km s<sup>-1</sup>. This covers the range of values this distribution is likely to have. This the same approach as adopted in our previous work \[Belczynski and Bulik, 1999\].

The binaries receive kicks for two reasons. First, the envelope of the supernova is lost from the system and it carries away some momentum. Thus even in the case when there is no kick velocity a binary achieves an additional velocity \[Blaauw, 1960\]. Second, if the supernova explosion is asymmetric both the newly formed compact object may receive a kick velocity which affects the orbit of the binary after the explosion as well as its center of mass velocity. The fate of a binary system in a supernova explosion depends on the value and direction of the kick velocity, on the orbital phase at which the explosion occurs, and on the parameters of the binary: the masses and orbital parameters $`a`$, and $`e`$.

We present the population of compact object binaries in the plane spanned by the center of mass velocity ofter the second supernova expolsion and time to merge in Figure 1. The orbital \[Blaauw, 1960\] effects are isolated and shown in the top left panel of Figure 1, where we present the results of the simulation with $`\sigma _v=0`$. There is a tail of long lived systems with lifetimes much longer than the Hubble time and small velocities, stretching outside of the boundaries of the plot to the lifetimes even of $`10^{20}`$years. These systems originally had large orbital separations, and hardly interacted in the course of their binary lifetime. In the case when there are no kicks the center of mass velocity of the comapct object binary depends on the amount of mass lost in the supernova explosion. In the extreme case of large mass loss, the center of mass velocity approaches the orbital velocity at the moment of supernova explosion, and it can never exceed it. The velocity of the system increases with increasing mass loss, however the systems that loose too much mass become unbound. This is why the lower part of the plot below $`t_{merge}10^8`$years is empty. With increasing the kick velocity also the typical velocity of a system increases and there appear short lived systems in tight orbits. They can now survive a large mass loss when the kick velocity has a favorable direction. Thus as the kick velocity is increased only the tightly bound systems (with short merger time) can survive the supernova explosions. Another effect of the kick velocity is that the long lived systems with $`t_{merge}`$ much longer than the Hubble time, which were present in the case $`\sigma _v=0`$km s<sup>-1</sup> disappear. The typical velocity of a system increases with the kick velocity. However, the population of the comapct merger binaries is not much affected when the kick velocity becomes large, e.g. changing the kick velocity distribution width from $`\sigma _v=0`$km s<sup>-1</sup> to $`200`$km s<sup>-1</sup> produces a much stronger effect than going from $`\sigma _v=400`$km s<sup>-1</sup> to $`800`$km s<sup>-1</sup>. Most of the systems are disrupted by such high velocities, and the surviving ones are only those for which the kick are not so large and have a favorable direction.

Another effect of increasing the kick velocity is that the typical lifetime of a system becomes smaller. When the kick velocity is large only very tight, and/or highly eccentric systems survive, hence the typical lifetime of compact object binaries decreases. It should be noted the typical center of mass velocity of the compact object binaries increases roughly linearly with the kick velocity, while the lifetime decreases approximately exponentially.

In Figure 1 we also plot following \[Bloom et al., 1998a\] the lines corresponding to the Hubble time (the dashed line), and we mark the region with the stars that will escape from a galactic potential. In order to escape a binary must satisfy the following conditions: (i) it has to have a velocity larger than the escape velocity, (ii) the distance $`vt_{merge}`$must be larger than the size of the galzxy. We also draw the line at $`t_{mrg}=15`$Myrs, to denote the systems that merger within the Hubble time. All the systems to the right of the solid line in Figure refvt have velocities above $`200`$km s<sup>-1</sup>, and live long enough to travel further than $`30`$kpc. We should also note that although each panel in Figure 1 contans $`10^3`$ systems, the production rate of compact object binaries decreases exponentially with the increasing kick velocity (see eq. 13 in Belczyński and Bulik \[Belczynski and Bulik, 1999\]).

In Figures 2 and 3 we present the cumulative distributions of the projected distance from the center of the host galaxy in case (i) and from their the place of birth in case (ii), respectively, of the systems that merge within the Hubble time. When the binaries propagate in the potential of a large galaxy the kick velocity only weakly influences the the distribution of the mergers. Below the radius of $`10`$kpc the distribution is determined by the potential well. This is where all short lived and slow systems merge. There exists however a tail of high velocity, long lived systems (see Figure 1) that manage to escape. The escaping fraction is a weak function of the kick velocity. Typically the number of systems that merge further than $`30`$kpc from the center of the host galaxy is $`30`$%, except for the unphysical case of no kick velocities when it drops below $`20`$%.

In the other extreme case of small host galaxies for which we neglect the gravitational potential the escaping fraction can be even larger. The escaping fraction decreases from $`80`$% for the kick velocity $`\sigma _v=200`$km s<sup>-1</sup> to about $`50`$% for $`\sigma _v=800`$km s<sup>-1</sup>. The reason for such behavior is clearly seen from Figure 1. In the product of center of mass velocity and time to merge the dominant role is played by the fast decrease of the time to merge with the increasing kick velocity $`\sigma _v`$.

These quantitative results are visualized in Figures 4 and 5. Here we show the distribution of $`10^3`$ mergers around massive galaxy and in the empty space. We are showing four panels that cover the scales from $`10`$kpc to $`10`$Mpc. In the case of the propagation in a massive galaxy we are showing the projection in the plane of the galactic disk so the effects of the rotational velocity and the asymmetry of the potential well are visible. Both calculations have been done for the case of the kick velocity distribution width $`\sigma _v=200`$km s<sup>-1</sup>.

## 4 Conclusions

We find that a significant fraction i.e. more than $`20`$ percent of compact object mergers take place outside of the host galaxies. The figures obtained for our case of no gravitational potential should be considered as an upper limit only. In contrast the observations show that the GRB afterglows lie on the host galaxies \[Hogg and Fruchter, 1998\]. However, our sample of the observed GRBs with afterglows is yet limited to the bursts longer than $`6`$s as BeppoSAX triggers on this timescale. Long bursts could be connected with the hypernovae-like events and therefore they are closely associated with the galaxies. Compact object mergers may be connected with the short bursts, although it has been argued that mass transfer in the coalescence of compact object may last much longer \[Portegies Zwart, 1998\].

Our results are consistent with the calculation by \[Bloom et al., 1998a\], for the case of a massive galaxy and the kick velocity distribution width $`\sigma _v=200`$km s<sup>-1</sup>. We have verified the dependence of the distribution of compact object mergers on $`\sigma _v`$. Our results show that for the case of massive galaxies the escaping fraction weakly depends on the distribution of kick velocities. We include also binaries with objects that are higher mass than the canonical $`1.4M_{}`$. These binaries are formed through accretion from a giant companion onto the neutron star. In this calculation the highest mass of a compact object is below $`2.5M_{}`$. The distribution of these more massive binaries is slightly more concentrated around the galaxies.

It has to be noted that there is a numnber of potential selection effects which may affect the results of this study. Assuming that compact objcet mergers are reponsible for GRBs there may be qualitative differences between the NS-NS mergers and NS-BH mergers. As indicated above, their spatial distribution around host galaxies is different. Also the typical timescale of the bursts may be different between these two classes. Gamma-ray bursts form two separate classes (long vs. short) with different brightness distributions and spectra (short burst are harder than the long ones). So far the study of afterglows has been possible only for the long bursts. It may be the case that compact object mergers are connected with the short bursts for which so far no information about the host galaxies exist.

The host galaxies have been identified in a long observational procedure: a gamma-ray burst lead to identification of a fading X-ray source, and then to discovery of the optical afterglow. Precise observations of the optical afterglows lead to the disory of host galaxies. There are bursts for which the X-ray or optical afterglows were not found. Since afterglows are usually connected with external shocks, gamma-ray bursts that take place outside of galaxies have much weaker afterglows because of the low density of intergalactic matter. Begelmann etal. Begelman and Meszaros \[Begelman et al., 1993\] argue that the afterglow emission depends scales only with the square root of the density of the outside medium. The mean external densities measured from the analysis of the known afterglow lightcurves are typically $`n0.03`$ cm<sup>-3</sup> \[Galama and Wijers, 1998\], while the intergalactic medium may be as rarified as $`10^6`$ cm<sup>-3</sup>. Hence the afterglow of a burst taking place outside a galaxy may be up to two orders of magnitude weaker than the one in a galaxy. It shows that there may be a strong preference against identification of host galaxies for the bursts that take place outside of galaxies.

Acknowledgments. This work has been funded by the following KBN grants: 2P03D01616 2P03D00911, 2P03D00415 and 2P03D01113, and also made use of the NASA Astrophysics Data System. TB is grateful for the hospitality of Ecole Polytechnique where this work was finished.