Evidence of a Past Merger of the Galactic Center Black Hole

We explore the evolution of SMBH mass and spin starting from three different types of BH seeds: those originating from Population III stars, gravitational runaways in star clusters, and direct collapse. Population III stars, forming within dark matter minihalos of $10^{5}$-$10^{6}$ $M_{\odot}$ from primordial gas are thought to be extremely massive, ranging from 10 to 1000 solar masses [44]. Their deaths lead to stellar mass BH seeds with masses ranging from 10-100 $M_{\odot}$[45, 46]. Alternatively, the formation of intermediate-mass black holes (IMBHs), ranging from $10^{3}$ to $10^{4}$ $M_{\odot}$, is hypothesized to occur through stellar mergers within the dense cores of stellar clusters in metal-poor protogalaxies[47, 48, 49, 50, 51, 52, 53, 54, 55], with each cluster having a mass of approximately $10^{5}$ $M_{\odot}$ [56, 57]. This process is facilitated by the dense gas environment, which increases the likelihood of stellar collisions. Lastly, some theories posit that the massive BH seeds can be formed from the direct collapse of clouds, resulting in massive BHs with mass ranges from $10^{4}$ to $10^{6}$ $M_{\odot}$[58, 59, 60]. Extended Data Figure 1 illustrates the mass distribution of the three types of BH seeds. For BHs from Population III stars, we employ stellar evolution with a Salpeter initial mass function. For gravitational runaways, we reference results from Ref. [57]. In the case of direct collapse, we assume a log-normal distribution centered at $\log_{10}(M_{\rm BH}/M_{\odot})=5$ with a standard deviation of 0.25.

Black hole accretion dynamics are generally classified into two types: cold flows[61, 62] , characterized by cooler, optically thick material at higher accretion rates, and hot flows[63, 64], which are hotter, optically thin, and occur at lower accretion rates. Hot flows, distinct for their lower radiative efficiency compared to the classical thin disk model, are often associated with the generation of astrophysical jets[65, 66, 67]. Such phenomena are observed in environments like low-luminosity active galactic nuclei and dormant black holes, including Sgr A*, the supermassive black hole at the center of the Milky Way.

Among various jet generation theories, the role of magnetic fields and black hole rotation is considered crucial. The Blandford-Znajek (BZ) model[20], in particular, has gained prominence. This model suggests that the rotational energy of the black hole itself powers the jets. The BZ model’s explanation aligns closely with observations and numerical simulations[68, 69, 70], indicating that relativistic jets primarily originate from the BZ process.

In the BZ model, a large-scale poloidal magnetic field penetrates the ergosphere of a rotating black hole and intersects the event horizon. The rotation of the black hole causes frame-dragging, which in turn generates a toroidal magnetic field and a Poynting flux. The essence of the BZ process is the possibility of an inward electromagnetic energy flux within the ergosphere that appears negative when measured from infinity. This negative flux crosses the event horizon, reducing the black hole’s mass-energy and angular momentum. Consequently, an outgoing jet emerges, Poynting-dominated, carrying away positive energy and angular momentum. The power of the jet in the BZ model can be described as[20],

$\displaystyle p_{\rm BZ}$

$\displaystyle\sim$

$\displaystyle 2.5\left(\frac{r_{g}}{r_{\bullet}}\right)^{2}a_{\bullet}^{2}\eta_{\rm BZ}^{2}\dot{M}c^{2},$

(1)

where $r_{g}$ and $r_{\bullet}$ represent the radii of Schwarzschild and Kerr black hole, respectivley, $a_{\bullet}=\frac{\mathbf{J}_{\bullet}c}{GM_{\bullet}^{2}}$ is the normalized black hole spin, and $\eta_{\rm BZ}=\Phi/\Phi_{\rm MAD}$ represents the ratio of the actual magnetic flux threading the black hole’s horizon ($\Phi$) to the maximum saturated magnetic flux ($\Phi_{\rm MAD}$). A value of $\eta_{\rm BZ}=1$ signifies that the magnetic flux from the magnetized accretion flow has saturated the black hole, leading to a magnetically arrested disk (MAD)[71, 72]. Conversely, a lower value of $\eta_{\rm BZ}$ suggests that the accretion disk is in a state of Standard And Normal Evolution (SANE)[73]. It is worth noting that Equation 1 applies for low to moderate spins. There are higher-order corrections as spin approaches 1[74, 75, 76]. The angular momentum carried out by the BZ jet then can be estimated by

$\displaystyle\dot{\mathbf{J}}_{\rm BZ}$

$\displaystyle=$

$\displaystyle p_{\rm BZ}/\Omega_{f}\hat{a_{\bullet}},$

(2)

where $\hat{a_{\bullet}}$ denotes the unit vector along the BH spin $a_{\bullet}$ vector, and

$\displaystyle\Omega_{f}$

$\displaystyle=$

$\displaystyle\frac{a_{\bullet}c}{2r_{\bullet}}=\frac{a_{\bullet}}{1+\sqrt{1-a_{\bullet}^{2}}}\frac{c}{2r_{g}}\,.$

(3)

Incorporated with the standard accretion model that excludes jet contributions, the equations governing the mass and angular momentum evolution of an accreting black hole can be expressed as follows[21]:

	$\displaystyle\frac{dM_{\bullet}}{dt}=\dot{M}E_{\rm ms}-p_{\rm BZ}/c^{2},$		(4)
	$\displaystyle\frac{d\mathbf{J_{\bullet}}}{dt}=\dot{M}J_{\rm ms}\hat{\mathbf{J}}_{d}-\dot{\mathbf{J}}_{\rm BZ},$		(5)

where $\hat{\mathbf{J}}_{d}$ is the unit vector along the disk angular momentum, $E_{\rm ms}$ and $J_{\rm ms}$ are the specific energy and specific angular momentum at the innermost radius $r_{\rm ISCO}$, which are defined as[77]

	$\displaystyle E_{\rm ms}$	$\displaystyle=$	$\displaystyle\frac{4\sqrt{r_{\rm ISCO}/r_{g}}-3a_{\bullet}}{\sqrt{3}r_{\rm ISCO}/r_{g}},$		(6)
	$\displaystyle J_{\rm ms}$	$\displaystyle=$	$\displaystyle\frac{GM_{\bullet}}{c}\frac{6\sqrt{r_{\rm ISCO}/r_{g}}-4a_{\bullet}}{\sqrt{3}\sqrt{r_{\rm ISCO}/r_{g}}}.$		(7)

For each accretion episode, Equations 4 and 5 are numerically solved using given parameterized values for $\dot{M}$ and $\eta_{\rm BZ}$.

During each accretion episode, we adopt a constant accretion rate $\dot{M}=\lambda_{\rm Edd}\dot{M}_{\rm Edd}$, where $\lambda_{\rm Edd}$ signifies the accretion efficiency, and $\dot{M}_{\rm Edd}=\frac{L_{\rm Edd}}{\epsilon c^{2}}$ represents the Eddington accretion rate. Here, the radiation efficiency $\epsilon$ is set to be 0.1. Accretion stops when the material within the self-gravitating radius $R_{\rm sg}$ is exhausted. The accretion time for each episode, $\tau_{\rm acc}$, can be approximated by $\frac{M_{d}}{\dot{M}}$, with $M_{d}$ being the total disk mass inside the self-gravitating radius $R_{\rm sg}$. The accretion time can be parameterized as[78],

$$\tau_{\rm acc}\sim\frac{M_{d}}{\dot{M}}=1.12\times 10^{6}\left(\frac{\alpha}{0.01}\right)^{-2/27}\left(\frac{\epsilon}{0.1}\right)^{22/27}\left(\frac{L}{0.1L_{\rm Edd}}\right)^{-22/27}M_{8}^{-4/27}{\rm yr}.$$

(8)

After an accretion episode concludes, the SMBH enters a dormant phase until the onset of the next episode. The interim period, known as the waiting time, is governed by the active galactic nucleus (AGN) duty cycle fraction $\tau_{\rm acc}/(\tau_{\rm acc}+\tau_{\rm wait})$, which is empirically observed to range from 1% to 10%.

The Lense-Thirring precession can result in the counter-alignment of a black hole’s spin vector $J_{\bullet}$ with the angular momentum vector $J_{d}$ of its accretion disc[10]. This counter-alignment process is contingent upon the condition that the angle $\theta$ between the vectors exceeds $\pi/2$, and that the angular momentum of the disc is less than twice that of the black hole, as $J_{d}<2J_{\bullet}$.

According to the analysis, if the angular momentum of the disc is less than twice that of the black hole, $J_{d}<2J_{\bullet}$, then for randomly oriented angular momentum vectors, the proportion $f_{c}$ of systems that will exhibit counter-alignment is given by the relationship:

$$f_{c}=\frac{1}{2}\left(1-\frac{J_{d}}{2J_{\bullet}}\right)$$

(9)

In scenarios where the disc’s angular momentum is less than that of the black hole $J_{d}<J_{\bullet}$, and accretion occurs through randomly oriented events, the black hole is likely to experience alternating episodes of spinup and spindown, a process referred to as chaotic accretion[10, 78, 79]. Conversely, in situations where $J_{d}>J_{\bullet}$, the black hole is expected to consistently gain angular momentum, or spin up, which is indicative of coherent accretion. This behavior would persist even if the black hole and the disc were initially counter-aligned.

In our accretion models, we characterize the extent of chaotic accretion using the Mises distribution function for the disk angular momentum vector $\mathbf{J}_{d}$[80]

$$p(\theta,k)=\frac{e^{k\cos\theta}}{2\pi J_{0}(k)},$$

(10)

where $J_{0}(k)$ is the zero-th order cylindrical Bessel function for the first kind. When $k=0$, the model yields a perfectly isotropic chaotic accretion. In contrast, a large $k>30$ results in nearly coherent accretion, where $\theta$ is consistently close to zero.

In the realm of binary black hole systems, those with black hole spins misaligned from the orbital angular momentum present significant challenges for analytical or semi-analytical modeling. The complex interactions between the spins, the orbital angular momentum, and each other lead to the precession of the system around the total angular momentum vector. To predict the final spin of the merged black hole, several faster approximate models have been developed[81, 82], covering waveforms[83] and remnant properties of merging black holes[84]. Those surrogate models for precessing binary black hole with generic spins and unequal masses are based on large sample of numerical relativity simulations, which are proven to be accurate in predicting the final spin magnitude and orientation with relative error at the order of $10^{-2}$ with 90th percentiles confidence[82]. With those surrogate models, we explore the probability that the SgrA* is the merger product of two SMBHs. We have adopted the latest NRSur7dq4EmriRemnant model[82], which is based on a seven-dimensional parameter space characterizing generically precessing binary black hole systems. This model incorporates the binary black hole mass ratio $q$, along with the two pre-merger spins $\mathbf{a}_{\bullet 1}$ and $\mathbf{a}_{\bullet 2}$.

In addition to the seven parameters described by the NRSur7dq4EmriRemnant model, our SMBH merger model includes an additional parameter, $\phi$, which represents the inclination angle of the SMBH binary angular momentum with respect to the LOS, as illustrated in Fig. 3. We explore 9 different inclinations ranging from $0^{\circ}$ to $180^{\circ}$, with 20^∘ equally spaced intervals. The spin magnitudes and orientations of $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ (with respect to the SMBH binary angular momentum) are randomly distributed, with magnitudes ranging from 0 to 1 to cover all possible pre-merger origins.

Recent studies of stellar surveys show specific chemical patterns and a retrograde motion in the nearby halo, suggesting a major past accretion event in the Galaxy. It has been found that the inner halo mainly consists of debris from a system, named Gaia-Enceladus[22], once larger than the Small Magellanic Cloud. Given the estimated mass ratio of 4:1, the merger between the Milky Way with Gaia-Enceladus around 10 billion years ago likely caused significant dynamical heating of the Galaxy’s early thick disk. Based on this observation, we fix the mass ratio of the merging black holes to be 4:1.

Owing to the 4:1 mass ratio, the final galactic plane is predominantly determined by the galactic plane of the pre-merger massive galaxy, which is our Milky Way. We assume that the final galactic plane aligns with the initial galactic plane of the more massive galaxy. The final spin in our merging model is computed using the surrogate model, which takes the previously described input parameters.

The cosmic population of SMBHBs is derived from the data of the Millennium-II N-body simulation [27], which has the resolution to distinguish dark matter halos down to $\sim 10^{8}M_{\odot}$. This dark matter framework is augmented with semianalytic models of galaxy evolution to trace the history of galaxies down to $\sim 10^{6}M_{\odot}$, which host central black holes as small as $\sim 10^{4}M_{\odot}$. All post-processing employs the same cosmological parameters as Millennium-II: $h=0.73$, $\Omega_{\Lambda}=0.75$, $\Omega_{M}=0.25$, and $\Omega_{k}=\Omega_{r}=0$. While these cosmological parameters are not the most current, variations of a few percent in $h$, $\Omega_{\Lambda}$, or $\Omega_{M}$ do not significantly affect our results.

We identify all mergers between galaxies, each containing an SMBH, totaling $N_{\bullet}=169435$ events throughout the Millennium-II simulation’s span. This was done by assuming an MBH occupation fraction of unity down to 10${}^{8}M_{\odot}$ stellar mass halos. This assumption slightly overestimates the SMBH binary merger rate at the lower mass end, as the MBH occupation fraction in these halos can drop to approximately $\sim$40%[85] for 10${}^{8}M_{\odot}$ stellar mass halos. We construct a 3D grid by logarithmically binning $\log_{10}(M_{1}/M_{\odot})$ and the mass ratio $q$ over the intervals $M_{1}\in[10^{5},10^{10}]M_{\odot}$ and $q\in[10^{-4},1]$, and linearly binning the redshift $z$ from 0 to 5 in line with the redshifts in Millennium-II snapshots. The redshift of each merger is matched with the snapshot where the progenitor galaxy is detected, assuming that the dynamical friction timescale for the SMBHs to form a close binary is shorter than the average 300 Myr between simulation snapshots.

By dividing the number of SMBHBs in each bin by the co-moving volume $V_{c}$ of the Millennium-II simulation, we approximate the differential SMBHB formation rate as

$$\frac{d^{4}N_{\bullet}}{dz\,dq\,dM_{1}\,dV_{c}}\,.$$

For those SMBHs that become a bound binary, the SMBHs start a long journey that bring them from separations of tens of kpc down to milli-parsec scales, below which they merge through GWs. The transition into the GW domain depends on energy exchanges with low-angular-momentum stars in the galaxy’s nucleus[86], or on the removal of angular momentum via gravitational torques from a gas disk, which can diminish their separation [87]. It might also involve a mix of these two mechanisms. Recent advancements in direct N-body simulations, Monte Carlo methods, and scattering experiments have provided a more hopeful outlook on what has been deemed the primary bottleneck in binary evolution for nearly four decades: the “final parsec problem”. This issue, described as the depletion of low-angular-momentum stars[86], seems less problematic. Findings suggest that the evolution of SMBHBs through stellar scattering potentially lead to a merger within less than approximately 1 billion years[88, 30, 89], especially when factors such as rotation, the triaxial shape of galaxies, and the granularity of stellar distribution are considered. SMBHB hardening through stellar scattering from pairing to the final coalescence can be described by[30]

	$\displaystyle\frac{da_{\rm b}}{dt}$	$\displaystyle=$	$\displaystyle\frac{da_{\rm b}}{dt}\bigg{|}_{3b}+\frac{da_{\rm b}}{dt}\bigg{|}_{\rm GW}$
		$\displaystyle=$	$\displaystyle-\frac{GH\rho_{\rm inf}}{\sigma_{\rm inf}}a^{2}_{\rm b}-\frac{64G^{3}M_{1}M_{2}(M_{1}+M_{2})f(e)}{5a_{\rm b}^{3}c^{5}},$

with

$$f(e)=\frac{1+(73/24)e^{2}+(37/96)e^{4}}{(1-e^{2})^{-7/2}},$$

(12)

where the first term represents the rate of binary hardening due to stellar scattering[90, 91], characterized by the dimensionless hardening rate $H$, stellar density $\rho_{\rm inf}$ and velocity dispersion $\sigma_{\rm inf}$ at the binary influence radius $r_{\rm inf}$. The influence radius is defined as the radius within which the enclosed stellar mass is twice the mass of the primary SMBH. The second term captures the effects of GW radiation[92, 93], which is characterized by the masses of the primary and secondary SMBHs in the binary ($M_{1}$ and $M_{2}$, respectively) and the binary’s eccentricity $e$. These terms together describe the evolution of the separation between the SMBHs under the combined influences of dynamical interactions with stars and energy loss due to GW radiation. The SMBHBs spend most of their lifetimes at the transition separation where $da_{b}/dt|_{3b}=da_{b}/dt|_{\rm GW}$, which is

$$a_{\rm GW}=\left(\frac{64G^{2}\sigma_{\rm inf}M_{1}M_{2}(M_{1}+M_{2})f(e)}{5c^{5}H\rho_{\rm inf}}\right)^{1/5}.$$

(13)

The corresponding timescale at this separation is

$$\tau_{\rm delay}\sim\frac{\sigma_{\rm inf}}{GH\rho_{\rm inf}a_{\rm GW}},$$

(14)

which can be used to estimate the time from SMBHB formation to final coalescence. The dimentionless hardening rate $H$ is estimated to be[30] $\sim 10$, and the velocity dispersion $\sigma_{\rm inf}$ can be estimated from the M-$\sigma$ relation [24]

$\displaystyle M_{9}$

$\displaystyle=$

$\displaystyle 0.309(\frac{\sigma_{\rm inf}}{200{\rm km\cdot s^{-1}}})^{4.38},$

(15)

where $M_{9}=M_{\bullet}/(10^{9}M_{\odot})$ and the stellar density at the influence radius can be obtained from the nuclei stellar cluster model[94]

$$\rho(r)=\frac{(3-\gamma)M_{*}}{4\pi}\frac{r_{0}}{r^{\gamma}(r+r_{0})^{4-\gamma}},$$

(16)

where $M_{*}\sim(1.84M_{9}^{0.86}\times 10^{11})M_{\odot}$ denotes the total mass of the stars in the bulge and $\gamma$ represents the inner cusp slope. The influence radius can be obtain by equations where[94, 95],

	$\displaystyle r_{\rm inf}$	$\displaystyle=$	$\displaystyle\frac{r_{0}}{(\frac{M_{*}}{2M_{\bullet}})^{1/(3-\gamma)}-1},$		(17)
	$\displaystyle r_{0}$	$\displaystyle=$	$\displaystyle\frac{4}{3}R_{\rm eff}(2^{1/(3-\gamma)}-1),$		(18)
	$\displaystyle R_{\rm eff}/{\rm pc}$	$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{ll}{\rm max}(2.95M_{*,6}^{0.596},34.8M_{*,6}^{0.399}),&{\quad\rm elliptical\quad galaxies},\\ 2.95M_{*,6}^{0.596},&{\quad\rm spiral\quad galaxies}.\end{array}\right.$		(21)

For each pair of SMBHs derived from the Millennium-II simulation, we can estimate the coalescence time starting from binary formation and calculate the redshift at which they ultimately merge using the following equation

	$\displaystyle\tau_{\rm delay}$	$\displaystyle=$	$\displaystyle\int_{z^{\prime}}^{z}\frac{dt}{d\bar{z}}d\bar{z}$		(22)
		$\displaystyle=$	$\displaystyle\frac{1}{H_{0}}\int_{z^{\prime}}^{z}\frac{d\bar{z}}{(1+\bar{z})\sqrt{\Omega_{r}(1+\bar{z})^{4}+\Omega_{M}(1+\bar{z})^{3}+\Omega_{k}(1+\bar{z})^{2}+\Omega_{\Lambda}}}\,,$

where $z$ denotes the redshift at which the SMBHB forms, and $z^{\prime}$ represents the redshift at which the SMBHB coalesces. If the coalescence delay time, $\tau_{\rm delay}$, is sufficiently long, it could result in an unphysical, negative $z^{\prime}$. This indicates that the SMBHB is still in the hardening phase. Such shrinking/stalled SMBHBs are excluded from our dataset. Then the differential SMBHB merger rate can be obtained as

$$\mathcal{R}_{\rm merger}(z^{\prime},q,M_{1})=\frac{d^{4}N_{\rm merger}}{dV_{c}\,dq\,dM_{1}\,dt}\bigg{|}_{z^{\prime}>0}\,.$$

(23)

The emission of gravitational waves during the merger of two SMBHs imparts a recoil, or ’kick’, to the resulting SMBH[96, 97]. If Sgr A* underwent a major merger in the past, the recoil velocity from this merger would cause Sgr A* to follow an elliptical orbit. After receiving this kick, the movement of the SMBH undergoes three distinct phases[39]. Initially, the SMBH oscillates back and forth with decreasing amplitude, losing energy through dynamical friction with each pass through the galactic core. This motion can be precisely modeled using dynamical friction. As the magnitude of motion decreases to about the core radius size, the SMBH and the core begin to oscillate around their common center of gravity, with these oscillations diminishing more slowly than the initial dynamical friction phase. In the final phase, the SMBH stabilizes and reaches thermal equilibrium with the surrounding stars. For the Milky Way, the stellar velocity dispersion at the Galactic center is approximately 100 km/s, leading to an equilibrium state of Brownian motion for Sgr A* at about 0.2 km/s[98]. Observations of Sgr A*’s peculiar motion perpendicular to the galactic plane[36, 37, 38] are consistent with this prediction, suggesting that Sgr A* is currently positioned quietly at the center of the Milky Way. However, the peculiar motion in the galactic plane is measured to be above the Brownian motion velocity. Therefore, the existence of a smaller BH or IMBH[99, 100, 101, 102, 103, 104, 105, 106, 107, 108] has not been completely ruled out. The most recent study posits a strong mass constraint on the secondary BH, indicating it to be less than 2000 $M_{\odot}$[109].

We employ the NRSur7dq4Remnant model to calculate the kick velocity of the final SMBH in our merger scenarios, as illustrated in the left panel of Extended Data Figure 8. Given that a 4:1 major merger with a spin parameter greater than or equal to 0.7 can accurately reproduce the spin of Sgr A*, it is likely that Sgr A* would receive a substantial kick velocity, ranging from hundreds to thousands of kilometers per second. By interpolating the results of Ref. [39](model A1 in their figure 15), we calculate the time required for Sgr A* to settle into the Brownian motion phase, with peculiar motion velocity $\sim$ 0.2 km/s. As demonstrated in the right panel of Extended Data Figure 8, for the entire parameter regime we explored, the recoil velocity Sgr A* acquired from the major merger could be sufficiently dissipated by the stars within the Milky Way’s core over the 9 billion years following the merger, consistent with the observations.