Future merger of the Milky Way with the Andromeda galaxy and the fate of their supermassive black holes
Riccardo Schiavi, Roberto Capuzzo-Dolcetta, Manuel Arca-Sedda, Mario Spera
AAstronomy & Astrophysics manuscript no. main © ESO 2021February 23, 2021
Future merger of the Milky Way with the Andromeda galaxy and thefate of their supermassive black holes
Riccardo Schiavi , Roberto Capuzzo-Dolcetta , Manuel Arca-Sedda , and Mario Spera , , Dipartimento di Fisica, Università degli Studi di Roma “La Sapienza”, P.le Aldo Moro, 2, I-00185 Rome, Italye-mail: [email protected] Astronomisches Rechen Intitut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstrasse 12-14, D-69120 Heidelberg,Germany Università di Padova, Dipartimento di Fisica e Astronomia, I-35131 Padova, Italy INFN, Sezione di Padova, I-35131 Padova, Italy Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Northwestern University, Evanston, IL 60208, USA
ABSTRACT
Our Galaxy and the nearby Andromeda galaxy (M31) are the most massive members of the Local Group, and they seem to be a boundpair, despite the uncertainties on the relative motion of the two galaxies. A number of studies have shown that the two galaxies willlikely undergo a close approach in the next 4 − N -body simulations to model this interaction to shed light on thefuture of the Milky Way - Andromeda system and for the first time explore the fate of the two supermassive black holes (SMBHs) thatare located at their centers. We investigated how the uncertainties on the relative motion of the two galaxies, linked with the initialvelocities and the density of the di ff use environment in which they move, a ff ect the estimate of the time they need to merge and form“Milkomeda”. After the galaxy merger, we follow the evolution of their two SMBHs up to their close pairing and fusion. Upon thefiducial set of parameters, we find that Milky Way and Andromeda will have their closest approach in the next 4.3 Gyr and merge overa span of 10 Gyr. Although the time of the first encounter is consistent with other predictions, we find that the merger occurs later thanpreviously estimated. We also show that the two SMBHs will spiral in the inner region of Milkomeda and coalesce in less than 16.6Myr after the merger of the two galaxies. Finally, we evaluate the gravitational-wave emission caused by the inspiral of the SMBHs,and we discuss the detectability of similar SMBH mergers in the nearby Universe ( z ≤
2) through next-generation gravitational-wavedetectors.
Key words. galaxies: interaction – galaxies: supermassive black holes – galaxies: Local Group – galaxies: kinematics and dynamics– gravitational waves – methods: numerical
1. Introduction
The Milky Way (MW) and the Andromeda galaxy (M31) are thetwo main members of the Local Group, which contains morethan 80 galaxies and has a total mass of roughly 3 − × M (cid:12) (González et al. 2014; van der Marel et al. 2012a). The futureevolution of the Local Group is essentially driven by the dynam-ics of our Galaxy and M31, and it can be considered a promis-ing study object to investigate the processes of galaxy formationand evolution. Although the physical and dynamical propertiesof the MW-M31 system are rather uncertain, it is likely that theLocal Group is gravitationally bound and decoupled from thegeneral cosmic expansion, and also that the two galaxies willnot escape the collision and the final merger. However, the timeat which this merger will occur is still a matter of debate. Themain purpose of this work and of our previous studies (Schiaviet 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 ap-proach will likely occur in < ff erent initial conditions in all the cited works. Usinga more recent estimation of the proper motion of M31, van derMarel et al. (2019) have obtained a time for the first approachequal to 4 . ff ected by a ratherhigh level of uncertainty. This is due mainly to the presence ofextended dark matter halos. We have chosen to adopt the valuesestimated in Klypin et al. (2002), defined as the virial masses atradius r , where the galactic density is 200 times the criticaldensity ρ ≈ . × − kg / m (according to the measurementsof the Hubble constant H by Huang et al. (2020) and Aghanimet al. (2018)): M MW = . × M (cid:12) and M M = . × M (cid:12) .Another source of uncertainty is our poor knowledge of the ac-tual size of the two galactic halos. The radial extent of the halo inequilibrium models of galaxies developed by Kuijken & Dubin-ski (1995) is in the range of 21 −
73 disk scale lengths. The ratioof the dark halo virial radius and the galaxy e ff ective radius fall Article number, page 1 of 8 a r X i v : . [ a s t r o - ph . GA ] F e b & A proofs: manuscript no. main in the same range in the studies by Jiang et al. (2019), Somervilleet al. (2017), Huang et al. (2017), and Kravtsov (2013). How-ever, there is evidence that some of the gaseous circumgalacticmedium (CGM) extends to even larger distances (Shull 2014). Inother words, we do not know with su ffi cient precision where agalaxy actually ends, and in our specific case, whether the MilkyWay and Andromeda are already partially overlapping or if theyare still well separated. For this reason, as discussed in Section2, we have decided to set the halo cuto ff radii of the two galaxiesat the respective tidal radii. Moreover, in Section 5.1 we demon-strate that the time of the merger does not depend on galactichalos that are more extended than 80 disk scale lengths. Be-cause it is evident that the edge of each galaxy gradually fades inthe IGM, we cannot ignore the e ff ect of this di ff use medium instudying the interaction. The density of the IGM is known to befour to six times the critical density ρ (Chamaraux & Tadokoro1971; Cox & Loeb 2008), but by performing several simulations,we obtained that even a small variation in this parameter coulda ff ect the merger time substantially.Our knowledge of the proper motion of M31 relative to us ismainly obtained through redshift measurement. This gained usan accurate estimate of the only radial component of the rela-tive velocity vector of M31: V r ≈
120 km / s (Binney & Tremaine1987). The tangential component has been inferred by studyingthe motion of the satellite galaxies of M31 (Loeb et al. 2005;van der Marel & Guhathakurta 2008) or by the Hubble SpaceTelescope (HST) and GAIA observations of sources behind M31(van der Marel et al. 2012b, 2014, 2016). These estimates spanfrom a minimum of V t ≈
17 km / s (van der Marel et al. 2012a) toa maximum of V t ≈
164 km / s (Salomon et al. 2016). The mostrecent estimate is that by van der Marel et al. (2019), who haveused GAIA DR2 to obtain V t = + − km / s. We referred to thisultimate measurement to fix the orientation of the relative ve-locity vector and its radial component ( V r = − . / s). Wediscuss the relation between the initial tangential velocity andthe time of the merger in Section 5.1.During the interaction at large scales, we are interested in fol-lowing the motion of the two SMBHs in the centers of the twogalaxies. It is well known that a compact object of mass M • MW = . × M (cid:12) (Gillessen et al. 2009), called SgrA*, is placed atthe center of the Milky Way. Even though the nucleus of M31seems to have a double or triple structure, there is a high proba-bility that it might host an SMBH of mass M • M = . × M (cid:12) (Bender et al. 2005). After the merger of the two host galaxies,their SMBHs are expected to form a binary that will shrink overtime through gravitational encounters with field stars. We firstexplore the future evolution of some of the nearest SMBHs andtheir eventual coalescence in the nucleus of the galactic mergerremnant. In Section 5.2 we discuss the time required for theSMBHs to merge and the amount of energy radiated throughgravitational waves (GWs).One of the most e ff ective ways to model galaxy interactions isthe integration of the N -body problem. While tree codes can sim-ulate a large number of collisionless particles in galaxy mergersimulations, a collisional direct summation N -body code withfewer particles is required in this study to follow the SMBH dy-namics. Direct N -body codes are highly reliable but computa-tionally expensive: this clearly places a limit on the number ofparticles involved in the simulations, and prevents us from re-solving large and small scales at the same time with good ac-curacy. For this reason, as we discuss in Section 3, we chose tosplit the whole study into two parts: in the first, we simulate thegalaxy interaction at large scales, and in the second, we focus onthe analysis of the orbital decay of the SMBH binary.
2. Galactic model and initial conditions
Our galaxies were modeled by combining three di ff erent com-ponents: an exponential disk, a spherical bulge, and a halo, thelatter two with a Hernquist (Hernquist 1990) density profile. Wecombined these three components with the command magalie in the NEMO code (Teuben 1995), which guarantees the stabilityof the whole system.The two galaxies have the structure presented in Klypin et al.(2002) and in Widrow & Dubinski (2005), that is, nearly thesame as was used by Cox & Loeb (2008). The main structureparameters are summarized in Table 1.
Milky Way AndromedaScale radius of the disk ( kpc ) . . Core radius of the bulge . . Core radius of the halo . . Cuto ff radius of the halo . . Mass of the disk ( M (cid:12) ) . × . × Mass of the bulge . . Mass of the halo . . Total mass ( M (cid:12) ) . × . × Table 1.
Values of characteristic parameters used in our simulations.Where it is not specified, the lengths are in units of the scale radius ofthe disk R d and masses are in units of the mass of the disk M d . The cuto ff radii of the halos were chosen as the tidal radii ofthe two galaxies, computed in the point-mass approximation: M MW M M = (cid:32) r t , MW R − r t , MW (cid:33) and r t , M = R − r t , MW , (1)where r t , MW and r t , M are the two tidal radii and R is the currentseparation between the two galaxies. A snapshot of our galacticmodel is shown in Fig. 1.We placed an SMBH in the center of each galaxy, which in thefirst part of our study was modeled as a particle with a mass of0.001 times the mass of the whole galaxy. This ratio implies thatthe mass of our SMBHs is significantly higher than the observedmasses, but this is not relevant for the dynamics of the two galax-ies until merging because the two SMBHs are essentially passiveguests of the hosting galaxies at this phase. We therefore madethis choice because it was the best setting allowed by our numer-ical resolution. We used a number of particles not greater than N = . × , and this constrains the mass of the single parti-cle. An SMBH mass of one thousandth of the mass of the galaxyis therefore a good compromise between the properties of thegalaxies and the comparison with an ordinary particle. However,during the simulation of the collision at large scales, the two par-ticles that represent the SMBHs only have the purpose of betteridentifying the two galactic centers and of observing their rela-tive distance at the end of the merger process.The Milky Way and Andromeda start to interact at the cur-rent distance of 780 kpc (McConnachie et al. 2005; Ribaset al. 2005), and their spin vectors are oriented at (0 ◦ ; − ◦ )and (240 ◦ ; − ◦ ) in Galactic coordinates, respectively (Dubin-ski et al. 1996; Raychaudury & Lynden-Bell 1989). To betterdisplay the dynamics of the galaxy binary system, we chose areference frame where the x-y plane coincides with the plane ofthe motion. The initial configuration of the two galaxies is shownin Fig. 2. Article number, page 2 of 8chiavi et al.: Future merger of the MW and M31 and the fate of their SMBHs -400-200 0 200 400-400 -200 0 200 400 y ( k p c ) x (kpc) -50-25 0 25 50-50 -25 0 25 50 y ( k p c ) x (kpc) Fig. 1.
Our model of the Milky Way. The three components (disk, bulge,and halo) are shown in di ff erent colors. The lower panel is a zoom intothe innermost region. 𝑦𝑧 𝑥 𝑀𝑊 𝑺 !" 𝑺 ! 𝑀31𝑽𝑽 % 𝑽 & Fig. 2.
Diagram of the initial configuration of the two galaxies in thechosen reference frame.
3. Methods
Owing to the computational complexity of the simulations, wedecided to divide the study into two parts. The first examines thegalaxy interaction at large scales and has the main purpose ofdetermining the true time required to form Milkomeda and itsfinal density profile. In this section we use a number of particlesof N = . × for the simulations with the fiducial set of pa-rameters and N = . × for all the others. The numericalintegration of the N-body equations of motion was implementedwith the HiGPUs code (Capuzzo Dolcetta et al. 2013). This pro- gram is based on a sixth-order Hermite integration scheme withblock time steps and directly computes the mutual force betweeneach pair of particles, exploiting a parallelization of CPUs andGPUs. Owing to the high performance of the
HiGPUs code, werepeated the simulation several times, changing two parametersthat were linked with the initial conditions and the external envi-ronment: the tangential component of the initial relative velocity V t , and the density of the IGM ρ . This allowed us to investigatethe correlation between the time of the merger and the galacticdynamical properties, together with the e ff ect of the dynamicalfriction exerted by the surrounding di ff use medium.In the second part we further study the evolution of the SMBHbinary that formed after the galaxy merger. Taking the simula-tion with the highest resolution and the fiducial set of parame-ters, we obtained the Milkomeda density profile and the velocitydispersion and modeled the galactic center as an analytic distri-bution of matter around the SMBH binary. To simulate the evolu-tion of the binary, we used a modified version of the ARWV code(Mikkola & Tanikawa 1999; Mikkola & Merritt 2008), whichintegrates the equations of motion taking in account the e ff ectof the dynamical friction exerted by a di ff use background dur-ing the first phases of the orbital evolution, the post-Newtonian(PN) terms when the binary shrinks enough to reach the GWemission regime, and the e ff ect of the spins of the merging ob-jects (Arca-Sedda & Capuzzo-Dolcetta 2019; Chassonnery et al.2019). To connect the new simulation to the previous one, theSMBHs orbit was reproduced with the same geometry as foundat the galaxy merger, while the masses of the two objects wereset to those known from observational evidences. Through thismethod, we can study the orbital decay of our binary down tosmall spatial scales and infer the coalescence time, the evolutionof the semimajor axis, the eccentricity, and the power emitted inthe form of GW.
4. Intergalactic medium and dynamical friction
The presence of the IGM a ff ects the time of the galaxy interac-tion through the extraction of orbital energy and angular momen-tum. We used HiGPUs to simulate the dynamics of the galaxycollision in di ff erent environments. To take the e ff ect of the IGMinto account, we modified the HiGPUs code by adding a dy-namical friction term in the equation of motion of each particleaccording to the Chandrasekhar formula (Binney & Tremaine1987), d r i dt = N (cid:88) j (cid:44) i Gm (cid:16) r j − r i (cid:17)(cid:18) (cid:15) + (cid:12)(cid:12)(cid:12) r j − r i (cid:12)(cid:12)(cid:12) (cid:19) / − π G ρ M ln Λ V c (cid:34) erf ( X ) − X √ π e − X (cid:35) V c , (2)with X = V c / √ σ , where σ is the IGM velocity dispersion.As usual for N -body codes, we introduced the softening param-eter (cid:15) to avoid the divergence of the Newton term at small dis-tances: it was fixed at (cid:15) =
500 pc for ordinary particles and (cid:15) =
50 pc for the two black holes. In the Chandrasekhar termwe considered ρ as the density of the IGM, and M and V c as themass and velocity of the galactic core that the particle belongs to.Unlike the classical Chandrasekhar formula, which describes thee ff ect of the dynamical friction on each star, depending on thestellar mass and velocity, in our case, each particle, which has Article number, page 3 of 8 & A proofs: manuscript no. main the same mass m , feels the same frictional force as all the oth-ers. This force changes with time during the galaxy interaction,but at any moment, is the same for every particle. This choice issuggested by the need of describing the collective e ff ect of thefriction on the motion of each galaxy as a whole.In all our simulations we fixed the Coulomb logarithm at ln Λ = σ = . / s,obtained from the equipartition of energy for a di ff use mediumat a temperature of T = × K (Cox & Loeb 2008).We compared the case with no IGM with three cases with di ff er-ent values of ρ : 1 . × − kg / m , the same as the critical density,4 . × − kg / m , which is the value estimated by Chamaraux &Tadokoro (1971), and 1 . × − kg / m , about 10 times the crit-ical density. As we show in Section 5.1, the time of the mergersignificantly changes for di ff erent IGM densities, especially forhigh initial velocities.
5. Results
For all initial conditions, the merger remnant Milkomeda re-sembles a giant elliptical galaxy with a density profile similarto those of the original two galaxies, as is shown in Fig. 3 for V t =
57 km / s and ρ = . × − kg / m . -2 den s i t y ( M / k p c ) r (kpc) Milky WayM31Milkomeda ⨀ Fig. 3.
Density profile of Milkomeda at time t = . V t =
57 km / s and ρ = . × − kg / m , and those of the MilkyWay and M31 at t = We obtain the time of the merger from the time evolutionof the distance between the centers of mass. The time of themerger is defined here as the time at which the separation is0.5% of its initial value.Before fixing the outermost edge of our galaxies at the respec-tive tidal radii, we investigated the dependence of the time ofthe merger T m on the cuto ff radius R h of the galactic halos. Inthe top panel of Fig. 4 we show that as expected, T m is stronglydependent on the galaxy extension only for low values of R h .For R h > R d , the timing of the process is no longer sensitiveto this parameter: from this value on, the two halos cover theentire distance between the two galaxies.We also found that when V t increases, T m rapidly increases, asis 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 accuraterelation between T m and V t : our best fit is T m ∝ V t .
15 20 25 30 35 0 50 100 150 200 T m ( G y r) R h /R d V t = 50 km/s T m ( G y r) V t (km/s) R h /R d = 70 Fig. 4.
Top: Time of the merger as a function of the cuto ff radius of thehalo R h , given in units of the scale radius of the disk R d and assumedequal for the two galaxies. The initial tangential velocity here is fixedat 50 km / s. Bottom: Correlation between the time of the merger and theinitial tangential velocity in the case of no IGM. The cuto ff radii of thehalos here are fixed at 70 R d . Repeating the simulation for di ff erent values of the IGM den-sity, we obtained that the presence of a di ff use medium speedsthe galaxy interaction up, especially in the case of large V t . InFig. 5 we plot the evolution of the distance between the two cen-ters of mass for four di ff erent values of ρ in the case of V t = / s (top panel) and the dependence of the time of the merger onthe IGM density for three di ff erent values of V t (bottom panel).Even though the time of the merger can significantly changewhen V t or ρ are varied, we note that the time of the first ap-proach is almost constrained in the interval 4 − ffi cient. The time of the first approachis very close to that obtained in the case of a pure radial fall oftwo point masses starting at a distance of 780 kpc with a relativeradial velocity of -115.7 km / s. After the first encounter, the IGMdensity instead plays a relevant role in the time for the comple-tion of the merger. This is mainly due to the enhanced speed ofthe two galaxies at the pericenter.Among the ensemble of the performed simulations, we refer tothe simulation with V t =
57 km / s and ρ = . × − kg / m as a fiducial case. According to our analysis, the Milky Way andAndromeda will reach their first approach in 4.3 Gyr and willmerge in 10.0 Gyr. The time for the first pericenter is close to Article number, page 4 of 8chiavi et al.: Future merger of the MW and M31 and the fate of their SMBHs the ∼ . ∼ . ∼ . ff erent to therecently measured velocity. Moreover, they used an IGM densitythat is slightly greater than we considered in our fiducial model. d i s t an c e ( k p c ) time (Gyr) V t = 57 km/s no IGM ρ = 10 -26 kg/m ρ = 4 × -26 kg/m ρ = 10 -25 kg/m T m ( G y r) ρ (10 -26 kg/m ) V t = 92 km/s V t = 57 km/s V t = 26 km/s Fig. 5.
Top: Separation between the two galaxies as a function of timefor three di ff erent values of IGM density, in the case of V t =
57 km / s.Bottom: Time of the merger as a function of the IGM density for threevalues of V t . We found that the distance between the two SMBHs located atthe galactic centers evolves in time, as was previously shownfor the two centers of mass. The only di ff erence is that after thegalaxy merger, the SMBH binary stalls at the same fixed dis-tance, independently of their initial velocity, as is shown in theFig. 6. This confirms the idea that the orbital decay of the bi-nary in the first phase simply follows the dynamics of the twostellar systems, but it later depends on the gravitational encoun-ters between the binary and the stars orbiting close to the galac-tic center. When the volume around the binary is depleted andthe binary orbit contains a total mass equal to or lower than thecombined mass of the two SMBHs, the binary stalls. Therefore,as expected, this interesting behavior is a function of the num-ber of particles that is involved in the simulation: the greater the number of particles, the smaller the stalling distance. Nev-ertheless, we note that as the numerical resolution increases at N > . × , the stalling distance reaches an asymptotic valueof ∼
100 pc, that is, about twice the softening parameter of thetwo SMBHs. Fig. 7 shows this behavior: for three di ff erent val-ues of N > . × , the stalling distance does not change.This might be the signal that we have reached the lower limit ofthe stalling distance that is allowed by our computational power.The density of stars around the binary rapidly drops to zero be-cause of the low numerical resolution, and this makes the en-ergy loss by dynamical friction ine ffi cient. However, because themain purpose of this first simulation is to reproduce the galaxyinteraction, we cannot expect to simultaneously correctly resolvethe dynamics at small scales. The stall shown in Fig. 7 is there-fore a direct e ff ect of the low spatial resolution and an indirecte ff ect of the sampling. d i s t an c e ( k p c ) time (Gyr) V t = 92 km/s V t = 57 km/s V t = 26 km/s Fig. 6.
Distance between the two SMBHs as a function of time for threedi ff erent initial tangential velocities. d i s t an c e ( k p c ) time (Gyr) N = 6.5 × N = 1.3 × N = 2.6 × Fig. 7.
Distance between the two SMBHs as a function of time for threedi ff erent numerical resolutions of the simulation. To follow the orbital decay of the SMBH binary, we used an-other code that improves the treatment of dynamical friction forthe second part of our study. We chose the simulation with thelargest number of particles ( N = . × ) and the fiducial set Article number, page 5 of 8 & A proofs: manuscript no. main of parameters ( V t =
57 km / s, ρ = . × − kg / m ) and cal-culated the density profile of Milkomeda at time t = . r h ≈ . γ = . . × M (cid:12) , to fit theinnermost region of the density profile of Milkomeda and re-produce the galactic center as an analytic external potential inthe ARWV code. The velocity dispersion in the innermost 500parsecs ( σ = . / s) was obtained from the outputs of HiGPUs code.The values of the masses of the two SMBHs were set to theproper values M • MW = . × M (cid:12) , and M • M = . × M (cid:12) (Gillessen et al. 2009; Bender et al. 2005), keeping as initial con-ditions those coming from the last computed orbit of the SMBHpair. The orbital integration with the ARWV code was performedin a reference frame in which the x-y plane coincides with theinitial plane of the motion.After loosing orbital energy owing to the interaction with the en-vironment, the binary becomes hard when the semimajor axisreaches the value defined in Merritt et al. (2007), a h = q (1 + q ) r h , (3)where q = M • MW / M • M = .
03 is the SMBH mass ratio. In ourmodel, r h ≈ . a h ≈ . a h .Through the comparison with the case in which the calcula-tion does not take the PN terms into account, we can infer thatmost of the orbital decay is due to the dynamical friction. Wedetermined that the only interaction with the background starscan carry the binary at a distance of a few times the character-istic Schwarzschild radius of the binary, which in our case is R s = . × − pc. In the absence of PN terms, the binarywould start to slowly shrink over a time of tens of million ofyears after this. The emission of gravitational waves rapidly re-duces the binary orbital energy and speeds up the decay: in afew thousand years, the distance plummets from few times R s to zero. According to our results, the merger between the twoSMBHs occurs 16.6 Myr after the formation of Milkomeda, inthe 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 theeccentricity e and the inverse semi-major axis 1 / a , respectively.The eccentricity starts at a value of 0.7 and drops to zero, as thebinary shrinks and circularizes.As expected, when the binary becomes hard, the hardening rate,defined by Merritt et al. (2007) as s = ddt (cid:32) a (cid:33) , (4)suddenly increases from ∼ − Myr − to ∼ × pc − Myr − in the last phases before the merger.However, because the external environment in ARWV is not mod-eled with an ensemble of particles, we cannot reproduce the or-bital energy loss due to the stochastic gravitational encounters of -5 -4 -3 -2 -1
0 2 4 6 8 10 12 14 16 1810 -15 -10 -5 d i s t an c e ( p c ) e m i tt ed po w e r ( L ) time (Myr)rel.dist. a GW.power -5 -4 -3 -2 -1
0 2 4 6 8 10 12 14 16 1810 -15 -10 -5 d i s t an c e ( p c ) e m i tt ed po w e r ( L ) time (Myr)rel.dist. a GW.power 10 -5 -4 -3 -2 -1
0 2 4 6 8 10 12 14 16 1810 -15 -10 -5 d i s t an c e ( p c ) e m i tt ed po w e r ( L ) time (Myr)rel.dist. a GW.power ⊙ Fig. 8.
Evolution of the relative distance between the two SMBHs (blueline) and the semimajor axis a of their orbit (red line), together withthe power emitted in the form of GWs (green line). The time is countedstarting from the galaxy merger. -4 -2
0 2 4 6 8 10 12 14 16 18 a - ( p c - ) , e time (Myr) a -1 e -5 -4 -3 -2 -1
0 2 4 6 8 10 12 14 16 1810 -15 -10 -5 d i s t an c e ( p c ) e m i tt ed po w e r ( L ) time (Myr)rel.dist. a GW.power
Fig. 9.
Inverse of the semimajor axis of the orbit of the SMBHs (in pc − )and its eccentricity as a function of time. The ordinate scale is the samefor a and e . 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 there-fore underestimated. We can obtain an estimate of the hardeningrate due to the energy exchange during the encounters, in the as-sumption of a background with a fixed and uniform density ρ anduniform velocity v , through the approximated formula (Quinlan1996; Gualandris et al. 2016) s = G ρ v H , (5)where H is a dimensionless hardening coe ffi cient with a nearlyconstant value of ∼
16 for hard binaries. In our case, this is s ≈ .
21 pc − Myr − , in accordance with those obtained in simi-lar simulations by Khan et al. (2018). 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 energyloss rate, averaged over one orbital period, according to the for-
Article number, page 6 of 8chiavi et al.: Future merger of the MW and M31 and the fate of their SMBHs mula (16) of Peters & Mathews (1963), (cid:104) P (cid:105) = G m m ( m + m ) c a (cid:0) − e (cid:1) / (cid:32) + e + e (cid:33) . (6)The emitted power progressively increases when the semi-major axis drops below a h and reaches a maximum of about10 L (cid:12) just before the merger.Through a numerical integration of the Peters formula, we ob-tained the amount of energy that is radiated away during the pro-cess. In Fig. 10 we compare the results of this integration withthe energy loss calculated by the ARWV code through the PN ap-proximation. The two curves agree well, with a fractionary log-arithmic variation below ∼ J.We repeated the simulation with the
ARWV code for two ex-
15 15.5 16 16.5 17 ene r g y ( J ) time (Myr)from eq. 8from ARWV
15 15.5 16 16.5 17 ene r g y ( J ) time (Myr)from eq. 8from ARWV Fig. 10.
Energy radiated during the last phases of the SMBH merger.The green line shows the mean emitted energy, obtained from the inte-gration of Eq. (6), and red points show the emitted energy as output ofthe
ARWV code. treme configurations of the SMBH spin vectors: one for paral-lel spins, and the other for antiparallel spins. As expected, thetime of the merger is not a ff ected by the spin orientation, but wefound that the final merger remnant gains a di ff erent recoil ve-locity. In the case of parallel spins, the recoil velocity of the rem-nant is v r = . / s, while for antiparallel spins, it is v r = . / s. The magnitude of the recoil velocity in both cases is small,mainly because of the low SMBH mass ratio (Healy & Lousto2018; Chassonnery & Capuzzo-Dolcetta 2020). The final SMBHis thus unable to escape from the center of the galaxy becausethe central escape velocity for our model is v e = . × km / s.Therefore, the giant elliptical galaxy Milkomeda will continueto host an SMBH in its center, as obtained also by Arca Seddaet al. (2019).We also investigated the possibility of observing a GW signalfrom a merger between similar SMBHs in the near Universe. Thefrequency-characteristic strain evolution is shown in Fig. 11 andoverlaps the sensitivity curve of di ff erent GW detectors, such asthe Pulsar Timing Array (PTA, Hobbs et al. (2010)), the SquareKilometer Array (SKA, Johnston et al. (2007)), the Laser Inter-ferometer Space Antenna (LISA, Amaro-Seoane et al. (2017)),the Deci-Hertz Interferometer Gravitational wave Observatory(DECIGO, Kawamura et al. (2011)), and the µ Ares microhertzdetector (Sesana et al. 2019). A comparison between our mod-eled signal and the detector sensitivities shows that mergers sim-ilar to the one we expect to witness in Milkomeda can be bright sources in ground-based detectors such as the PTA, or in the nextdecade, the SKA, provided that they take place roughly within1 Mpc. However, farther away in space, it becomes clear thatthe only possibility of observing this type of SMBH mergers isusing space-borne detectors. Mergers that occur up to redshift z (cid:46) µ Ares detector concept design. − − − − − − − Frequency [Hz]10 − − − − − − − − C h a r a c t e r i s t i c S tr a i n L I S A D E C I G O µ A r e s P T A S K A MilkomedaObservation time = 4 yr, a = 1980 AU - e = 0 . z = 0 , D L = 780 kpc z = 0 . z = 1 . z = 2 . Fig. 11.
Characteristic strain of a GW signal emitted by a similar SMBHmerging binary as a function of frequency for four di ff erent redshifts.The GW signal is computed starting from a semimajor axis of a = e = . t =
6. Discussion and conclusions
We studied the future evolution of the system composed of theMilky Way and M31 (Andromeda galaxy) in relation with thegravitational e ff ects due to the intergalactic background. Thiscan shed light not only on the forthcoming dynamics of our spe-cific galaxy system, but also in general on the correlation be-tween the galaxy interaction timing and the environmental prop-erties, the galactic structure, and the initial conditions. We usedthe high-precision N -body HiGPUs code, properly modified toaccount for the friction exerted on the bodies in the galaxies bythe di ff use background, to calculate the evolution of the galac-tic hosts and their central SMBHs until the merger of the twogalaxies. At this stage, we used the galaxy density, the veloc-ity distribution, and the SMBH orbital parameters for a furtherdynamical simulation. With the ARWV code in its most recentversion, which considers PN treatment and the BH spins, wereproduced the orbital decay of the SMBH binary in the post-merger galactic background. The aim was not simply to estimatethe likelihood and future time of an eventual merger, but also todetermine the fate of the two massive black holes hosted at thetwo galactic centers, whose estimated mass is 4 . × M (cid:12) inour Galaxy and 1 . × M (cid:12) in M31, giving a mass ratio 0 . − Article number, page 7 of 8 & A proofs: manuscript no. main for the background density, the dimension of the halos, andthe initial velocity.2. The time for the completion of the two galaxy merger in-creases significantly with the relative velocity transversecomponent, which is an ill-determined observational quan-tity. According to the most recent estimates, however we canconclude that the MW and M31 will merge in ∼
10 Gyr.3. The dimensions of the galactic halos play an important rolein the merger time: larger halos cause a significant orbitalenergy dissipation and accelerate the decay, at least until thedistance between the two galaxies is on the same order as thesize of the halos.4. As expected, due to the collisionless nature of the encounterand merger, the overall post-merger density profile is notvery di ff erent from a mere mass average of the two profilesof the MW and M31.5. After the merger, the two SMBHs were left orbiting on amutual orbit of eccentricity e ∼ . a ∼
160 pc, which stalled because the resolution of the N -bodysimulation is insu ffi cient.6. The following fate of the SMBH pair was followed by a PNsimulation that showed how e ffi ciently (in less than 17 Myr)dynamical friction braking leads the two SMBHs to the so-called hard binary phase, when subsequent orbital decay isgiven by energy dissipation by GW emission down to thefinal merger and recoil kick.7. When we also considered antiparallel spins for the SMBHs,the recoil kick velocity was below 25 km / s (two orders ofmagnitude lower than the central escape velocity), whichleaves the BH remnant confined in the inner potential wellof the galaxy.8. Types of SMBH orbital decays similar to those studied hereshow a very high power of GW emission because of the highmasses of the BH involved. This high power is shifted towardvery low frequency. This GW emission would in principlebe observable only with future GW ground-based detectorssuch as the PTA and the SKA or with space interferometerssuch as LISA, but the redshift range for the detection shouldbe 1 ≤ z ≤ References