EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 642 (October 2020) A&A, 642 (2020) A30 Full HTML
Free Access
Issue
A&A
Volume 642, October 2020
Article Number A30
Number of page(s) 8
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202038674
Published online 01 October 2020

© ESO 2020

1. Introduction

The Milky Way (MW) and the Andromeda galaxy (M 31) are the two main members of the Local Group, which contains more than 80 galaxies and has a total mass of roughly 3 − 5 × 1012 M⊙ (González et al. 2014; van der Marel et al. 2012a). The future evolution of the Local Group is essentially driven by the dynamics of our Galaxy and M 31, and it can be considered a promising study object to investigate the processes of galaxy formation and evolution. Although the physical and dynamical properties of the MW-M 31 system are rather uncertain, it is likely that the Local Group is gravitationally bound and decoupled from the general cosmic expansion, and also that the two galaxies will not escape the collision and the final merger. However, the time at which this merger will occur is still a matter of debate. The main purpose of this work and of our previous studies (Schiavi et al. 2019a,b) is to shed light on this topic.

According to some previous simulations (Dubinski et al. 1996; Cox & Loeb 2008; van der Marel et al. 2012b), the first close approach will likely occur in < 4 Gyr, even though the two galaxies have different initial conditions in all the cited works. Using a more recent estimation of the proper motion of M 31, van der Marel (2019) have obtained a time for the first approach equal to 4.5 Gyr. Almost the same result can also be obtained when the so-called “timing argument” is employed, which was introduced in the pioneering work by Kahn & Woltjer (1959). In the timing argument, MW and M 31 are considered as point masses on a radial orbit: they started their motion at the Big Bang, and after decoupling from the Hubble flow, began to approach one another. Even though the timing argument allowed obtaining an estimate of the total mass of the Local Group that is compatible with the estimate obtained with other methods (e.g., Klypin et al. 2002; Widrow & Dubinski 2005), it is unable to take the complexity of the dynamics of the galaxy interaction into account. The time needed for the completion of the whole merger process is highly sensitive not only to the masses of the two galaxies, but also to their proper motion and to the density of the intergalactic medium (IGM) in which they move. All the estimates of the mass of the two galaxies are affected by a rather high level of uncertainty. This is due mainly to the presence of extended dark matter halos. We have chosen to adopt the values estimated in Klypin et al. (2002), defined as the virial masses at radius r200, where the galactic density is 200 times the critical density ρ0 ≈ 1.0 × 10−26 kg m−3 (according to the measurements of the Hubble constant H0 by Huang et al. 2020; Planck Collaboration VI 2020): MMW = 1.0 × 1012 M⊙ and MM 31 = 1.6 × 1012 M⊙. Another source of uncertainty is our poor knowledge of the actual size of the two galactic halos. The radial extent of the halo in equilibrium models of galaxies developed by Kuijken & Dubinski (1995) is in the range of 21 − 73 disk scale lengths. The ratio of the dark halo virial radius and the galaxy effective radius fall in the same range in the studies by Jiang et al. (2019), Somerville et al. (2018), Huang et al. (2017), and Kravtsov (2013). However, there is evidence that some of the gaseous circumgalactic medium (CGM) extends to even larger distances (Shull 2014). In other words, we do not know with sufficient precision where a galaxy actually ends, and in our specific case, whether the Milky Way and Andromeda are already partially overlapping or if they are still well separated. For this reason, as discussed in Sect. 2, we have decided to set the halo cutoff radii of the two galaxies at the respective tidal radii. Moreover, in Sect. 5.1 we demonstrate that the time of the merger does not depend on galactic halos that are more extended than 80 disk scale lengths. Because it is evident that the edge of each galaxy gradually fades in the IGM, we cannot ignore the effect of this diffuse medium in studying the interaction. The density of the IGM is known to be four to six times the critical density ρ0 (Chamaraux & Tadokoro 1971; Cox & Loeb 2008), but by performing several simulations, we obtained that even a small variation in this parameter could affect the merger time substantially.

Our knowledge of the proper motion of M 31 relative to us is mainly obtained through redshift measurement. This gained us an accurate estimate of the only radial component of the relative velocity vector of M 31: Vr ≈ 120 km s−1 (Binney & Tremaine 1987). The tangential component has been inferred by studying the motion of the satellite galaxies of M 31 (Loeb et al. 2005; van der Marel & Guhathakurta 2008) or by the Hubble Space Telescope (HST) and Gaia observations of sources behind M 31 (van der Marel et al. 2012b; van der Marel & Kallivayalil 2014; van der Marel & Sahlmann 2016). These estimates span from a minimum of Vt ≈ 17 km s−1 (van der Marel et al. 2012a) to a maximum of Vt ≈ 164 km s−1 (Salomon et al. 2016). The most recent estimate is that by van der Marel (2019), who have used Gaia DR2 to obtain V t = 57 − 31 + 35 $ V_{\mathrm{t}}=57_{-31}^{+35} $ km s−1. We referred to this ultimate measurement to fix the orientation of the relative velocity vector and its radial component (Vr = −115.7 km s−1). We discuss the relation between the initial tangential velocity and the time of the merger in Sect. 5.1.

During the interaction at large scales, we are interested in following the motion of the two SMBHs in the centers of the two galaxies. It is well known that a compact object of mass M MW • = 4.31 × 10 6 M ⊙ $ M_{\mathrm{MW}}^{\bullet}=4.31\times10^6\,{M_{\odot}} $ (Gillessen et al. 2009), called SgrA*, is placed at the center of the Milky Way. Even though the nucleus of M 31 seems to have a double or triple structure, there is a high probability that it might host an SMBH of mass M M 31 • = 1.4 × 10 8 M ⊙ $ M_{\mathrm{M\,31}}^{\bullet}=1.4\times10^8\,{M_{\odot}} $ (Bender et al. 2005). After the merger of the two host galaxies, their SMBHs are expected to form a binary that will shrink over time through gravitational encounters with field stars. We first explore the future evolution of some of the nearest SMBHs and their eventual coalescence in the nucleus of the galactic merger remnant. In Sect. 5.2 we discuss the time required for the SMBHs to merge and the amount of energy radiated through gravitational waves (GWs).

One of the most effective ways to model galaxy interactions is the integration of the N-body problem. While tree codes can simulate a large number of collisionless particles in galaxy merger simulations, a collisional direct summation N-body code with fewer particles is required in this study to follow the SMBH dynamics. Direct N-body codes are highly reliable but computationally expensive: this clearly places a limit on the number of particles involved in the simulations, and prevents us from resolving large and small scales at the same time with good accuracy. For this reason, as we discuss in Sect. 3, we chose to split the whole study into two parts: in the first, we simulate the galaxy interaction at large scales, and in the second, we focus on the analysis of the orbital decay of the SMBH binary.

2. Galactic model and initial conditions

Our galaxies were modeled by combining three different components: an exponential disk, a spherical bulge, and a halo, the latter two with a Hernquist (Hernquist 1990) density profile. We combined these three components with the command magalie in the NEMO code (Teuben 1995), which guarantees the stability of the whole system.

The two galaxies have the structure presented in Klypin et al. (2002) and in Widrow & Dubinski (2005), that is, nearly the same as was used by Cox & Loeb (2008). The main structure parameters are summarized in Table 1.

Table 1.

Values of characteristic parameters used in our simulations.

The cutoff radii of the halos were chosen as the tidal radii of the two galaxies, computed in the point-mass approximation:

M MW M M 31 = ( r t , MW R − r t , MW ) 2 and r t , M 31 = R − r t , MW , $$ \begin{aligned} \frac{M_{\rm MW}}{M_{\rm M\,31}}=\left(\frac{r_{\rm t,MW}}{R-r_{\rm t,MW}}\right)^2\quad \mathrm{and}\quad r_{\rm t,M\,31}=R-r_{\rm t,MW}, \end{aligned} $$(1)

where rt, MW and rt, M 31 are the two tidal radii and R is the current separation between the two galaxies. A snapshot of our galactic model is shown in Fig. 1.

thumbnail Fig. 1.

Our model of the Milky Way. The three components (disk, bulge, and halo) are shown in different colors. Lower panel: zoom into the innermost region.

We placed an SMBH in the center of each galaxy, which in the first part of our study was modeled as a particle with a mass of 0.001 times the mass of the whole galaxy. This ratio implies that the mass of our SMBHs is significantly higher than the observed masses, but this is not relevant for the dynamics of the two galaxies until merging because the two SMBHs are essentially passive guests of the hosting galaxies at this phase. We therefore made this choice because it was the best setting allowed by our numerical resolution. We used a number of particles not greater than N = 2.6 × 105, and this constrains the mass of the single particle. An SMBH mass of one thousandth of the mass of the galaxy is therefore a good compromise between the properties of the galaxies and the comparison with an ordinary particle. However, during the simulation of the collision at large scales, the two particles that represent the SMBHs only have the purpose of better identifying the two galactic centers and of observing their relative distance at the end of the merger process.

The Milky Way and Andromeda start to interact at the current distance of 780 kpc (McConnachie et al. 2005; Ribas et al. 2005), and their spin vectors are oriented at (0° ; − 90° ) and (240° ; − 30° ) in Galactic coordinates, respectively (Dubinski et al. 1996; Raychaudury & Lynden-Bell 1989). To better display the dynamics of the galaxy binary system, we chose a reference frame where the x − y plane coincides with the plane of the motion. The initial configuration of the two galaxies is shown in Fig. 2.

thumbnail Fig. 2.

Diagram of the initial configuration of the two galaxies in the chosen reference frame.

3. Methods

Owing to the computational complexity of the simulations, we decided to divide the study into two parts. The first examines the galaxy interaction at large scales and has the main purpose of determining the true time required to form Milkomeda and its final density profile. In this section we use a number of particles of N = 2.6 × 105 for the simulations with the fiducial set of parameters and N = 6.5 × 104 for all the others. The numerical integration of the N-body equations of motion was implemented with the HiGPUs code (Capuzzo Dolcetta et al. 2013). This program is based on a sixth-order Hermite integration scheme with block time steps and directly computes the mutual force between each pair of particles, exploiting a parallelization of CPUs and GPUs. Owing to the high performance of the HiGPUs code, we repeated the simulation several times, changing two parameters that were linked with the initial conditions and the external environment: the tangential component of the initial relative velocity Vt, and the density of the IGM ρ. This allowed us to investigate the correlation between the time of the merger and the galactic dynamical properties, together with the effect of the dynamical friction exerted by the surrounding diffuse medium.

In the second part we further study the evolution of the SMBH binary that formed after the galaxy merger. Taking the simulation with the highest resolution and the fiducial set of parameters, we obtained the Milkomeda density profile and the velocity dispersion and modeled the galactic center as an analytic distribution of matter around the SMBH binary. To simulate the evolution of the binary, we used a modified version of the ARWV code (Mikkola & Tanikawa 1999; Mikkola & Merritt 2008), which integrates the equations of motion taking in account the effect of the dynamical friction exerted by a diffuse background during the first phases of the orbital evolution, the post-Newtonian (PN) terms when the binary shrinks enough to reach the GW emission regime, and the effect of the spins of the merging objects (Arca-Sedda & Capuzzo-Dolcetta 2019; Chassonnery et al. 2019). To connect the new simulation to the previous one, the SMBHs orbit was reproduced with the same geometry as found at the galaxy merger, while the masses of the two objects were set to those known from observational evidences. Through this method, we can study the orbital decay of our binary down to small spatial scales and infer the coalescence time, the evolution of the semimajor axis, the eccentricity, and the power emitted in the form of GW.

4. Intergalactic medium and dynamical friction

The presence of the IGM affects the time of the galaxy interaction through the extraction of orbital energy and angular momentum. We used HiGPUs to simulate the dynamics of the galaxy collision in different environments. To take the effect of the IGM into account, we modified the HiGPUs code by adding a dynamical friction term in the equation of motion of each particle according to the Chandrasekhar formula (Binney & Tremaine 1987),

d 2 r i d t 2 = ∑ j ≠ i N G m ( r j − r i ) ( ϵ 2 + | r j − r i | 2 ) 3 / 2 − 4 π G 2 ρ M ln Λ V c 3 [ erf ( X ) − 2 X π e − X 2 ] V c , $$ \begin{aligned} \frac{d^2{\boldsymbol{r}}_i}{\mathrm{d}t^2}=&\sum \limits _{j\ne i}^N \frac{Gm\left({\boldsymbol{r}}_j-{\boldsymbol{r}}_i\right)}{\left(\epsilon ^2+\left|{\boldsymbol{r}}_j-{\boldsymbol{r}}_i\right|^2\right)^{3/2}} \nonumber \\&-\frac{4\pi G^2\rho M\ln {\Lambda }}{V_{\rm c}^3}\left[\mathrm{erf}{(X)}-\frac{2X}{\sqrt{\pi }}e^{-X^2}\right]{\boldsymbol{V}}_{\rm c}, \end{aligned} $$(2)

with X = V c / 2 σ $ X=V_{\mathrm{c}}/\sqrt{2}\sigma $, where σ is the IGM velocity dispersion.

As usual for N-body codes, we introduced the softening parameter ϵ to avoid the divergence of the Newton term at small distances: it was fixed at ϵ = 500 pc for ordinary particles and ϵ = 50 pc for the two black holes. In the Chandrasekhar term we considered ρ as the density of the IGM, and M and Vc as the mass and velocity of the galactic core that the particle belongs to. Unlike the classical Chandrasekhar formula, which describes the effect of the dynamical friction on each star, depending on the stellar mass and velocity, in our case, each particle, which has the same mass m, feels the same frictional force as all the others. This force changes with time during the galaxy interaction, but at any moment, is the same for every particle. This choice is suggested by the need of describing the collective effect of the friction on the motion of each galaxy as a whole.

In all our simulations we fixed the Coulomb logarithm at lnΛ = 5 and the velocity dispersion of the medium at σ = 86.2 km s−1, obtained from the equipartition of energy for a diffuse medium at a temperature of T = 3 × 105 K (Cox & Loeb 2008).

We compared the case with no IGM with three cases with different values of ρ: 1.0 × 10−26 kg m−3, the same as the critical density, 4.0 × 10−26 kg m−3, which is the value estimated by Chamaraux & Tadokoro (1971), and 1.0 × 10−25 kg m−3, about 10 times the critical density. As we show in Sect. 5.1, the time of the merger significantly changes for different IGM densities, especially for high initial velocities.

5. Results

5.1. Galaxy merger

For all initial conditions, the merger remnant Milkomeda resembles a giant elliptical galaxy with a density profile similar to those of the original two galaxies, as is shown in Fig. 3 for Vt = 57 km s−1 and ρ = 4.0 × 10−26 kg m3.

thumbnail Fig. 3.

Density profile of Milkomeda at time t = 10.2 Gyr in the case with Vt = 57 km s−1 and ρ = 4.0 × 10−26 kg m3, and those of the Milky Way and M 31 at t = 0.

We obtain the time of the merger from the time evolution of the distance between the centers of mass. The time of the merger is defined here as the time at which the separation is 0.5% of its initial value.

Before fixing the outermost edge of our galaxies at the respective tidal radii, we investigated the dependence of the time of the merger Tm on the cutoff radius Rh of the galactic halos. In the top panel of Fig. 4 we show that as expected, Tm is strongly dependent on the galaxy extension only for low values of Rh. For Rh >  80Rd, the timing of the process is no longer sensitive to this parameter: from this value on, the two halos cover the entire distance between the two galaxies.

thumbnail Fig. 4.

Top: time of the merger as a function of the cutoff radius of the halo Rh, given in units of the scale radius of the disk Rd and assumed equal for the two galaxies. The initial tangential velocity here is fixed at 50 km s−1. Bottom: correlation between the time of the merger and the initial tangential velocity in the case of no IGM. The cutoff radii of the halos here are fixed at 70 Rd.

We also found that when Vt increases, Tm rapidly increases, as is shown in the lower panel of Fig. 4 in the case of no IGM. It is interesting to note that there seems to be a quite accurate relation between Tm and Vt: our best fit is T m ∝ V t 5 $ T_{\rm m} \propto {V_{\rm t}}^5 $.

Repeating the simulation for different values of the IGM density, we obtained that the presence of a diffuse medium speeds the galaxy interaction up, especially in the case of large Vt. In Fig. 5 we plot the evolution of the distance between the two centers of mass for four different values of ρ in the case of Vt = 57 km s−1 (top panel) and the dependence of the time of the merger on the IGM density for three different values of Vt (bottom panel).

thumbnail Fig. 5.

Top: separation between the two galaxies as a function of time for three different values of IGM density, in the case of Vt = 57 km s−1. Bottom: time of the merger as a function of the IGM density for three values of Vt.

Even though the time of the merger can significantly change when Vt or ρ are varied, we note that the time of the first approach is almost constrained in the interval 4−5 Gyr. This means that the first part of the galaxy motion is nearly Keplerian because the orbital energy dissipation due to the friction exerted by the IGM is still not very efficient. The time of the first approach is very close to that obtained in the case of a pure radial fall of two point masses starting at a distance of 780 kpc with a relative radial velocity of −115.7 km s−1. After the first encounter, the IGM density instead plays a relevant role in the time for the completion of the merger. This is mainly due to the enhanced speed of the two galaxies at the pericenter.

Among the ensemble of the performed simulations, we refer to the simulation with Vt = 57 km s−1 and ρ = 4.0 × 10−26 kg m3 as a fiducial case. According to our analysis, the Milky Way and Andromeda will reach their first approach in 4.3 Gyr and will merge in 10.0 Gyr. The time for the first pericenter is close to the ∼4.5 Gyr found by van der Marel (2019), but they did not report any value for the time of the final merger. Cox & Loeb (2008) obtained the first approach at ∼2.8 Gyr and the merger at ∼5.4 Gyr. However, we have to consider that they started the simulations of the MW-M 31 interaction 5 Gyr in the past and reached a current transverse velocity that is very different to the recently measured velocity. Moreover, they used an IGM density that is slightly greater than we considered in our fiducial model.

5.2. SMBH merger

We found that the distance between the two SMBHs located at the galactic centers evolves in time, as was previously shown for the two centers of mass. The only difference is that after the galaxy merger, the SMBH binary stalls at the same fixed distance, independently of their initial velocity, as is shown in the Fig. 6. This confirms the idea that the orbital decay of the binary in the first phase simply follows the dynamics of the two stellar systems, but it later depends on the gravitational encounters between the binary and the stars orbiting close to the galactic center. When the volume around the binary is depleted and the binary orbit contains a total mass equal to or lower than the combined mass of the two SMBHs, the binary stalls. Therefore, as expected, this interesting behavior is a function of the number of particles that is involved in the simulation: the greater the number of particles, the smaller the stalling distance. Nevertheless, we note that as the numerical resolution increases at N >  5.0 × 104, the stalling distance reaches an asymptotic value of ∼100 pc, that is, about twice the softening parameter of the two SMBHs. Figure 7 shows this behavior: for three different values of N >  5.0 × 104, the stalling distance does not change. This might be the signal that we have reached the lower limit of the stalling distance that is allowed by our computational power. The density of stars around the binary rapidly drops to zero because of the low numerical resolution, and this makes the energy loss by dynamical friction inefficient. However, because the main purpose of this first simulation is to reproduce the galaxy interaction, we cannot expect to simultaneously correctly resolve the dynamics at small scales. The stall shown in Fig. 7 is therefore a direct effect of the low spatial resolution and an indirect effect of the sampling.

thumbnail Fig. 6.

Distance between the two SMBHs as a function of time for three different initial tangential velocities.

thumbnail Fig. 7.

Distance between the two SMBHs as a function of time for three different numerical resolutions of the simulation.

5.2.1. Orbital decay of the SMBH binary

To follow the orbital decay of the SMBH binary, we used another code that improves the treatment of dynamical friction for the second part of our study. We chose the simulation with the largest number of particles (N = 2.6 × 105) and the fiducial set of parameters (Vt = 57 km s−1, ρ = 4.0 × 10−26 kg m3) and calculated the density profile of Milkomeda at time t = 10.2 Gyr, soon after the merger, when the system has stabilized and the SMBH binary has formed. This last point occurs when the two objects approximately reach the so-called influence radius, defined in Merritt et al. (2007) as the radius of a sphere around the two SMBHs that contains twice the sum of their masses. In our case, this corresponds to rh ≈ 156.5 pc.

We used a Dehnen power law (Dehnen 1993), with γ = 0.8, scale length 50 kpc, and total mass 9.3 × 1013 M⊙, to fit the innermost region of the density profile of Milkomeda and reproduce the galactic center as an analytic external potential in the ARWV code. The velocity dispersion in the innermost 500 parsecs (σ = 203.5 km s−1) was obtained from the outputs of HiGPUs code.

The values of the masses of the two SMBHs were set to the proper values M MW • = 4.31 × 10 6 M ⊙ $ M_{\mathrm{MW}}^{\bullet}=4.31\times10^6\,{M_{\odot}} $, and M M 31 • = 1.4 × 10 8 M ⊙ $ M_{\mathrm{M\,31}}^{\bullet}=1.4\times10^8\,{M_{\odot}} $ (Gillessen et al. 2009; Bender et al. 2005), keeping as initial conditions those coming from the last computed orbit of the SMBH pair. The orbital integration with the ARWV code was performed in a reference frame in which the x − y plane coincides with the initial plane of the motion.

After loosing orbital energy owing to the interaction with the environment, the binary becomes hard when the semimajor axis reaches the value defined in Merritt et al. (2007),

a h = q ( 1 + q ) 2 r h 4 , $$ \begin{aligned} a_{\rm h}=\frac{q}{(1+q)^2}\frac{r_{\rm h}}{4}, \end{aligned} $$(3)

where q= M MW • / M M   31 • =0.03 $ q=M_{\rm MW}^\bullet/M_{\rm M\,31}^\bullet=0.03 $ is the SMBH mass ratio. In our model, rh ≈ 156.6 pc, and therefore ah ≈ 1.1 pc. Figure 8 shows the distance between the two SMBHs and the semimajor axis of their orbit, together with the emitted GW power, as a function of time. In absolute values, the slopes of the three lines start to increase rapidly when the separation approximately reaches ah.

thumbnail Fig. 8.

Evolution of the relative distance between the two SMBHs (blue line) and the semimajor axis a of their orbit (red line), together with the power emitted in the form of GWs (green line). The time is counted starting from the galaxy merger.

Through the comparison with the case in which the calculation does not take the PN terms into account, we can infer that most of the orbital decay is due to the dynamical friction. We determined that the only interaction with the background stars can carry the binary at a distance of a few times the characteristic Schwarzschild radius of the binary, which in our case is Rs = 1.38 × 10−5 pc. In the absence of PN terms, the binary would start to slowly shrink over a time of tens of million of years after this. The emission of gravitational waves rapidly reduces the binary orbital energy and speeds up the decay: in a few thousand years, the distance plummets from few times Rs to zero. According to our results, the merger between the two SMBHs occurs 16.6 Myr after the formation of Milkomeda, in the same range of timescales as was found by Khan et al. (2016) for similar SMBH mergers.

The blue and red curves in Fig. 9 show the evolution of the eccentricity e and the inverse semi-major axis 1/a, respectively. The eccentricity starts at a value of 0.7 and drops to zero, as the binary shrinks and circularizes.

thumbnail Fig. 9.

Inverse of the semimajor axis of the orbit of the SMBHs (in pc−1) and its eccentricity as a function of time. The ordinate scale is the same for a and e.

As expected, when the binary becomes hard, the hardening rate, defined by Merritt et al. (2007) as

s = d d t ( 1 a ) , $$ \begin{aligned} s=\frac{d}{\mathrm{d}t}\left(\frac{1}{a}\right), \end{aligned} $$(4)

suddenly increases from ∼4 pc−1 Myr−1 to ∼4 × 104 pc−1 Myr−1 in the last phases before the merger.

However, because the external environment in ARWV is not modeled with an ensemble of particles, we cannot reproduce the orbital energy loss due to the stochastic gravitational encounters of the binary with the field stars. Our estimation of s, and of |de/dt| as well, in the phases soon after the binary becomes hard is therefore underestimated. We can obtain an estimate of the hardening rate due to the energy exchange during the encounters, in the assumption of a background with a fixed and uniform density ρ and uniform velocity v, through the approximated formula (Quinlan 1996; Gualandris et al. 2016)

s = G ρ v H , $$ \begin{aligned} s=\frac{G\rho }{{ v}}H, \end{aligned} $$(5)

where H is a dimensionless hardening coefficient with a nearly constant value of ∼16 for hard binaries. In our case, this is s ≈ 0.21 pc−1 Myr−1, in accordance with those obtained in similar simulations by Khan et al. (2018).

5.2.2. Gravitational wave emission

The power emitted in the form of GWs as a function of time, shown with the green line in Fig. 8, was obtained as the energy loss rate, averaged over one orbital period, according to the formula (16) of Peters & Mathews (1963),

⟨ P ⟩ = 32 5 G 4 m 1 2 m 2 2 ( m 1 + m 2 ) c 5 a 5 ( 1 − e 2 ) 7 / 2 ( 1 + 73 24 e 2 + 37 96 e 4 ) · $$ \begin{aligned} \langle P \rangle =\frac{32}{5}\frac{G^4m_1^2m_2^2\left(m_1+m_2\right)}{c^5a^5\left(1-e^2\right)^{7/2}}\left(1+\frac{73}{24}e^2+\frac{37}{96}e^4\right)\cdot \end{aligned} $$(6)

The emitted power progressively increases when the semimajor axis drops below ah and reaches a maximum of about 1019 L⊙ just before the merger.

Through a numerical integration of the Peters formula, we obtained the amount of energy that is radiated away during the process. In Fig. 10 we compare the results of this integration with the energy loss calculated by the ARWV code through the PN approximation. The two curves agree well, with a fractionary logarithmic variation below ∼13%. The amount of energy emitted in the last phases of process is on the order of 1043 J.

thumbnail Fig. 10.

Energy radiated during the last phases of the SMBH merger. The green line shows the mean emitted energy, obtained from the integration of Eq. (6), and red points show the emitted energy as output of the ARWV code.

We repeated the simulation with the ARWV code for two extreme configurations of the SMBH spin vectors: one for parallel spins, and the other for antiparallel spins. As expected, the time of the merger is not affected by the spin orientation, but we found that the final merger remnant gains a different recoil velocity. In the case of parallel spins, the recoil velocity of the remnant is vr = 2.2 km s−1, while for antiparallel spins, it is vr = 24.8 km s−1. The magnitude of the recoil velocity in both cases is small, mainly because of the low SMBH mass ratio (Healy & Lousto 2018; Chassonnery & Capuzzo-Dolcetta, in prep.). The final SMBH is thus unable to escape from the center of the galaxy because the central escape velocity for our model is ve = 3.7 × 103 km s−1. Therefore, the giant elliptical galaxy Milkomeda will continue to host an SMBH in its center, as obtained also by Arca Sedda et al. (2019).

We also investigated the possibility of observing a GW signal from a merger between similar SMBHs in the near Universe. The frequency-characteristic strain evolution is shown in Fig. 11 and overlaps the sensitivity curve of different GW detectors, such as the Pulsar Timing Array (PTA, Hobbs et al. 2010), the Square Kilometer Array (SKA, Johnston et al. 2007), the Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2017), the Deci-Hertz Interferometer Gravitational wave Observatory (DECIGO, Kawamura et al. 2011), and the μAres microhertz detector (Sesana et al. 2019). A comparison between our modeled signal and the detector sensitivities shows that mergers similar to the one we expect to witness in Milkomeda can be bright sources in ground-based detectors such as the PTA, or in the next decade, the SKA, provided that they take place roughly within 1 Mpc. However, farther away in space, it becomes clear that the only possibility of observing this type of SMBH mergers is using space-borne detectors. Mergers that occur up to redshift z ≲ 2 might be observable by LISA, if with an allegedly low signal-to-noise ratio, but they might shine bright to a micro-hertz observatory such as the μAres detector concept design.

thumbnail Fig. 11.

Characteristic strain of a GW signal emitted by a similar SMBH merging binary as a function of frequency for four different redshifts. The GW signal is computed starting from a semimajor axis of a0 = 1980 AU and an eccentricity of e0 = 0.064, corresponding to a time of t = 4 yr before the merger. The sensibility curves of the main GW detectors are also shown.

6. Discussion and conclusions

We studied the future evolution of the system composed of the Milky Way and M 31 (Andromeda galaxy) in relation with the gravitational effects due to the intergalactic background. This can shed light not only on the forthcoming dynamics of our specific galaxy system, but also in general on the correlation between the galaxy interaction timing and the environmental properties, the galactic structure, and the initial conditions. We used the high-precision N-body HiGPUs code, properly modified to account for the friction exerted on the bodies in the galaxies by the diffuse background, to calculate the evolution of the galactic hosts and their central SMBHs until the merger of the two galaxies. At this stage, we used the galaxy density, the velocity distribution, and the SMBH orbital parameters for a further dynamical simulation. With the ARWV code in its most recent version, which considers PN treatment and the BH spins, we reproduced the orbital decay of the SMBH binary in the post-merger galactic background. The aim was not simply to estimate the likelihood and future time of an eventual merger, but also to determine the fate of the two massive black holes hosted at the twos galactic centers, whose estimated mass is 4.31 × 106 M⊙ in our Galaxy and 1.4 × 108 M⊙ in M 31, giving a mass ratio 0.03. We summarize our results below.

  1. The time evolution of the MW and M 31 orbits is such that the first close approach of the two galaxies will occur in 4 − 5 Gyr, with a weak dependence on the characteristics assumed for the background density, the dimension of the halos, and the initial velocity.

  2. The time for the completion of the two galaxy merger increases significantly with the relative velocity transverse component, which is an ill-determined observational quantity. According to the most recent estimates, however we can conclude that the MW and M 31 will merge in ∼10 Gyr.

  3. The dimensions of the galactic halos play an important role in the merger time: larger halos cause a significant orbital energy dissipation and accelerate the decay, at least until the distance between the two galaxies is on the same order as the size of the halos.

  4. As expected, due to the collisionless nature of the encounter and merger, the overall post-merger density profile is not very different from a mere mass average of the two profiles of the MW and M 31.

  5. After the merger, the two SMBHs were left orbiting on a mutual orbit of eccentricity e ∼ 0.7 and semimajor axis a ∼ 160 pc, which stalled because the resolution of the N-body simulation is insufficient.

  6. The following fate of the SMBH pair was followed by a PN simulation that showed how efficiently (in less than 17 Myr) dynamical friction braking leads the two SMBHs to the so-called hard binary phase, when subsequent orbital decay is given by energy dissipation by GW emission down to the final merger and recoil kick.

  7. When we also considered antiparallel spins for the SMBHs, the recoil kick velocity was below 25 km s−1 (two orders of magnitude lower than the central escape velocity), which leaves the BH remnant confined in the inner potential well of the galaxy.

  8. Types of SMBH orbital decays similar to those studied here show a very high power of GW emission because of the high masses of the BH involved. This high power is shifted toward very low frequency. This GW emission would in principle be observable only with future GW ground-based detectors such as the PTA and the SKA or with space interferometers such as LISA, but the redshift range for the detection should be 1 ≤ z ≤ 2.

On the basis of our results, we are now able to determine the most feasible scenario of the future of our own Galaxy and its central SMBH. Our new estimate of the time required for the completion of the MW-M 31 merger means that the life of the Local Group is slightly longer than previously believed. A final result is that unfortunately, the Sun will not live long enough to witness the formation of Milkomeda and will therefore not be part of the new galaxy.

References

  1. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, ArXiv e-prints [arXiv:1702.00786] [Google Scholar]
  3. Arca-Sedda, M., & Capuzzo-Dolcetta, R. 2019, MNRAS, 483, 152 [NASA ADS] [CrossRef] [Google Scholar]
  4. Arca Sedda, M., Berczik, P., Capuzzo-Dolcetta, R., et al. 2019, MNRAS, 484, 520 [Google Scholar]
  5. Bender, R., Kormendy, J., Bower, G., et al. 2005, ApJ, 631, 280 [NASA ADS] [CrossRef] [Google Scholar]
  6. Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, NJ: Princeton University Press) [Google Scholar]
  7. Capuzzo Dolcetta, R., Spera, M., Punzo, D., et al. 2013, J. Comput. Phys., 236, 580 [NASA ADS] [CrossRef] [Google Scholar]
  8. Chamaraux, P., & Tadokoro, M. 1971, PASJ, 23, 117 [NASA ADS] [Google Scholar]
  9. Chassonnery, P., Capuzzo-Dolcetta, R., & Mikkola, S. 2019, ArXiv e-prints [arXiv:1910.05202] [Google Scholar]
  10. Cox, T. J., & Loeb, A. 2008, MNRAS, 386, 461 [NASA ADS] [CrossRef] [Google Scholar]
  11. Dehnen, W. 1993, MNRAS, 265, 250 [NASA ADS] [CrossRef] [Google Scholar]
  12. Dubinski, J., Mihos, J. C., Hernquist, L., et al. 1996, ApJ, 462, 576 [NASA ADS] [CrossRef] [Google Scholar]
  13. Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 1075 [NASA ADS] [CrossRef] [Google Scholar]
  14. González, R. E., Kravtsov, A. V., Gnedin, N. Y., et al. 2014, ApJ, 793, 91 [NASA ADS] [CrossRef] [Google Scholar]
  15. Gualandris, A., Read, J. I., Dehnen, W., et al. 2016, MNRAS, 464, 2301 [NASA ADS] [CrossRef] [Google Scholar]
  16. Healy, J., & Lousto, C. O. 2018, Phys. Rev. D, 97, 084002 [NASA ADS] [CrossRef] [Google Scholar]
  17. Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
  18. Hobbs, G., Archibald, A., Arzoumanian, Z., Backer, D., et al. 2010, Classical Quantum Gravity, 27, 084013 [NASA ADS] [CrossRef] [Google Scholar]
  19. Huang, K.-H., Fall, S. M., Ferguson, H. C., et al. 2017, ApJ, 838, 6 [NASA ADS] [CrossRef] [Google Scholar]
  20. Huang, C. D., Riess, A. G., Hoffmann, S. L., et al. 2020, ApJ, 889, 5 [NASA ADS] [CrossRef] [Google Scholar]
  21. Jiang, F., Dekel, A., Kneller, O., et al. 2019, MNRAS, 488, 4801 [CrossRef] [Google Scholar]
  22. Johnston, S., Bailes, M., Bartel, N., et al. 2007, PASA, 24, 174 [NASA ADS] [CrossRef] [Google Scholar]
  23. Kahn, F. D., & Woltjer, L. 1959, ApJ, 130, 705 [NASA ADS] [CrossRef] [Google Scholar]
  24. Kawamura, S., Ando, M., Seto, N., et al. 2011, Classical Quantum Gravity, 28, 094011 [NASA ADS] [CrossRef] [Google Scholar]
  25. Khan, F. M., Fiacconi, D., Mayer, L., Berczik, P., & Just, A. 2016, ApJ, 828, 73 [NASA ADS] [CrossRef] [Google Scholar]
  26. Khan, F. M., Capelo, P. R., Mayer, L., et al. 2018, ApJ, 868, 97 [Google Scholar]
  27. Klypin, A., Zhao, H., Somerville, R., et al. 2002, ApJ, 573, 597 [Google Scholar]
  28. Kravtsov, A. V. 2013, ApJ, 764, L31 [NASA ADS] [CrossRef] [Google Scholar]
  29. Kuijken, K., & Dubinski, J. 1995, MNRAS, 277, 1341 [CrossRef] [Google Scholar]
  30. Loeb, A., Reid, M. J., Brunthaler, A., et al. 2005, ApJ, 633, 894 [NASA ADS] [CrossRef] [Google Scholar]
  31. McConnachie, A. W., Irwin, M. J., Ferguson, A. M. N., et al. 2005, MNRAS, 356, 979 [NASA ADS] [CrossRef] [Google Scholar]
  32. Merritt, D., Mikkola, S., Szell, A., et al. 2007, ApJ, 671, 53 [NASA ADS] [CrossRef] [Google Scholar]
  33. Mikkola, S., & Merritt, D. 2008, ApJ, 135, 2398 [NASA ADS] [CrossRef] [Google Scholar]
  34. Mikkola, S., & Tanikawa, K. 1999, MNRAS, 310, 745 [CrossRef] [Google Scholar]
  35. Peters, P. C., & Mathews, J. 1963, Phys. Rev., 131, 435 [NASA ADS] [CrossRef] [Google Scholar]
  36. Quinlan, G. D. 1996, New Astron., 1, 35 [NASA ADS] [CrossRef] [Google Scholar]
  37. Raychaudury, S., & Lynden-Bell, D. 1989, MNRAS, 240, 195 [NASA ADS] [CrossRef] [Google Scholar]
  38. Ribas, I., Jordi, C., Vilardell, F., et al. 2005, ApJ, 635, L37 [NASA ADS] [CrossRef] [Google Scholar]
  39. Salomon, J.-B., Ibata, R. A., Famaey, B., et al. 2016, MNRAS, 456, 4432 [CrossRef] [Google Scholar]
  40. Schiavi, R., Capuzzo-Dolcetta, R., & Arca Sedda, M. 2019a, ArXiv e-prints [arXiv:1906.02982] [Google Scholar]
  41. Schiavi, R., Capuzzo-Dolcetta, R., Arca Sedda, M., & Spera, M. 2019b, Proc. Int. Astron. Union, 14, 161 [CrossRef] [Google Scholar]
  42. Sesana, A., Korsakova, N., Arca Sedda, M., et al. 2019, ArXiv e-prints [arXiv:1908.11391] [Google Scholar]
  43. Shull, J. M. 2014, ApJ, 784, 142 [NASA ADS] [CrossRef] [Google Scholar]
  44. Somerville, R. S., Behroozi, P., Pandya, V., et al. 2018, MNRAS, 473, 2714 [NASA ADS] [CrossRef] [Google Scholar]
  45. Teuben, P. 1995, ASP Conf. Ser., 77, 398 [NASA ADS] [Google Scholar]
  46. van der Marel, R. P. 2019, ApJ, 872, 24 [NASA ADS] [CrossRef] [Google Scholar]
  47. van der Marel, R. P., & Guhathakurta, P. 2008, ApJ, 678, 187 [NASA ADS] [CrossRef] [Google Scholar]
  48. van der Marel, R. P., & Kallivayalil, N. 2014, ApJ, 781, 121 [NASA ADS] [CrossRef] [Google Scholar]
  49. van der Marel, R. P., & Sahlmann, J. 2016, ApJ, 832, L23 [CrossRef] [Google Scholar]
  50. van der Marel, R. P., Besla, G., Cox, T. J., et al. 2012a, ApJ, 753, 9 [NASA ADS] [CrossRef] [Google Scholar]
  51. van der Marel, R. P., Fardal, M., Besla, G., et al. 2012b, ApJ, 753, 8 [NASA ADS] [CrossRef] [Google Scholar]
  52. Widrow, L. M., & Dubinski, J. 2005, ApJ, 631, 838 [Google Scholar]

All Tables

Table 1.

Values of characteristic parameters used in our simulations.

In the text

All Figures

thumbnail Fig. 1.

Our model of the Milky Way. The three components (disk, bulge, and halo) are shown in different colors. Lower panel: zoom into the innermost region.
In the text
thumbnail Fig. 2.

Diagram of the initial configuration of the two galaxies in the chosen reference frame.
In the text
thumbnail Fig. 3.

Density profile of Milkomeda at time t = 10.2 Gyr in the case with Vt = 57 km s−1 and ρ = 4.0 × 10−26 kg m3, and those of the Milky Way and M 31 at t = 0.
In the text
thumbnail Fig. 4.

Top: time of the merger as a function of the cutoff radius of the halo Rh, given in units of the scale radius of the disk Rd and assumed equal for the two galaxies. The initial tangential velocity here is fixed at 50 km s−1. Bottom: correlation between the time of the merger and the initial tangential velocity in the case of no IGM. The cutoff radii of the halos here are fixed at 70 Rd.
In the text
thumbnail Fig. 5.

Top: separation between the two galaxies as a function of time for three different values of IGM density, in the case of Vt = 57 km s−1. Bottom: time of the merger as a function of the IGM density for three values of Vt.
In the text
thumbnail Fig. 6.

Distance between the two SMBHs as a function of time for three different initial tangential velocities.
In the text
thumbnail Fig. 7.

Distance between the two SMBHs as a function of time for three different numerical resolutions of the simulation.
In the text
thumbnail Fig. 8.

Evolution of the relative distance between the two SMBHs (blue line) and the semimajor axis a of their orbit (red line), together with the power emitted in the form of GWs (green line). The time is counted starting from the galaxy merger.
In the text
thumbnail Fig. 9.

Inverse of the semimajor axis of the orbit of the SMBHs (in pc−1) and its eccentricity as a function of time. The ordinate scale is the same for a and e.
In the text
thumbnail Fig. 10.

Energy radiated during the last phases of the SMBH merger. The green line shows the mean emitted energy, obtained from the integration of Eq. (6), and red points show the emitted energy as output of the ARWV code.
In the text
thumbnail Fig. 11.

Characteristic strain of a GW signal emitted by a similar SMBH merging binary as a function of frequency for four different redshifts. The GW signal is computed starting from a semimajor axis of a0 = 1980 AU and an eccentricity of e0 = 0.064, corresponding to a time of t = 4 yr before the merger. The sensibility curves of the main GW detectors are also shown.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.