© The Authors 2024

1. Introduction

The formation of supermassive black holes (BHs) likely involves initial seeds undergoing accretion and merger episodes (Volonteri et al. 2021). Empirical relations between the BH mass, the stellar dispersion velocity, and the mass of the galactic bulge (McConnell & Ma 2013) suggest the occurrence of common evolutionary processes by which the growth of the BH is linked to that of its host galaxy and influences its properties. This interplay is most evident for active galactic nuclei (AGNs): material falling toward the bottom of the galaxy potential well forms an accretion disk around the supermassive BH that, in turn, shapes the surrounding environment.

Most galaxies are expected to undergo at least one major merger during their evolution (O’Leary et al. 2021), and this is especially true for the brightest cluster galaxies sitting at the center of clusters (Gerosa & Sesana 2015). Galactic mergers plays a fundamental role in the formation of cosmic structures and provide new influxes of gas, ultimately leading to an increased rate of star formation and the powering of AGNs. Following a galaxy merger, the supermassive BHs hosted by the two galaxies are expected to form binaries and eventually merge – a process that likely involves a variety of mechanisms, including dynamical friction, stellar scattering, gas-assisted migration, and gravitational wave (GW) emission (Begelman et al. 1980; Colpi 2014). A potential strategy for inferring the occurrence of supermassive BH mergers besides direct GW probes (Amaro-Seoane et al. 2017) is to detect post-merger electromagnetic signatures (e.g., Madau & Quataert 2004; Bonning et al. 2007; Blecha et al. 2016).

Anisotropic emission of GWs during the late BH inspiral and merger imparts a proper velocity, or “kick”, to the BH remnant (Bonnor & Rotenberg 1961; Peres 1962; Bekenstein 1973). The kick velocity is independent of the total mass of the binary and is typically 𝒪(100 km s⁻¹) (see, e.g., Gerosa et al. 2018). However, some fine-tuned spin configurations give rise to “superkicks” of 𝒪(100 km s⁻¹) (González et al. 2007; Campanelli et al. 2007). Such velocities are comparable to, if not higher than, the escape speed of the galactic host itself (Merritt et al. 2004).

There are two main observational strategies used to identify recoiling supermassive BHs; both assume that some of the nearby material remains bound to the BH after the merger (Loeb 2007; Volonteri & Madau 2008):

• Astrometric surveys are used to identify off-nuclear AGNs, targeting objects with a projected position that is somewhat displaced from the galactic center (e.g., Madau & Quataert 2004).

• Spectroscopic searches for recoiling BHs target the predicted Doppler shift between the so-called broad lines, which originate from material bound to the BH, and the narrow lines, which instead track the host galaxy (e.g., Bonning et al. 2007).

Even though these signatures are not unique, several recoiling supermassive BH candidates have been proposed (Komossa et al. 2008; Batcheldor et al. 2010; Robinson et al. 2010; Civano et al. 2012; Steinhardt et al. 2012; Comerford & Greene 2014; Koss et al. 2014; Lena et al. 2014; Markakis et al. 2015; Barrows et al. 2016; Shapovalova et al. 2016; Chiaberge et al. 2017; Kim et al. 2017, 2018; Kalfountzou et al. 2017; Pesce et al. 2018, 2021; Skipper & Browne 2018; Hogg et al. 2021), some of which have since been ruled out (Decarli et al. 2014; Li et al. 2024). The most promising candidate to date is arguably the quasar 3C 186 (Chiaberge et al. 2017, 2018). This source presents both an evident astrometric displacement of ∼10 kpc between the quasar and the isophotal center of the galaxy stellar distribution and a ∼2000 km s⁻¹ Doppler shift between broad and narrow lines. Using submillimeter interferometry, Castignani et al. (2022) reported the detection of a large reservoir of molecular gas (≳10¹⁰ M_⊙) orbiting the center of the observed galaxy; at the same time, cold gas was not detected at the location and redshift of the quasar, in line with the recoiling BH scenario¹.

Morishita et al. (2022) point out that the projected direction of the jet (which presumably tracks the spin of the BH) in 3C 186 is essentially perpendicular to the projected separation between the AGN and the photometric center of the galaxy (which presumably tracks the direction of the kick). If 3C 186 is a post-merger BH, its kick and spin appear perpendicular to each other. This is puzzling, at least at first sight, because general relativity strongly predicts that large kicks of 𝒪(1000) km s⁻¹ tend to be parallel to the BH spin (Gerosa et al. 2018). In this paper we show that this new piece of observational evidence together with detailed relativistic modeling allows for an unprecedented characterization of 3C 186. Compared to the previous study by Lousto et al. (2017), our analysis includes a complete modeling of all vector components of spin and kick, which are crucial to solving the geometry of the system. Our statistical methods are also more advanced, similar to those used by Paynter & Thrane (2023) for the recoiling BH candidate E1821+643 (Shapovalova et al. 2016). The key advantage of studying 3C 186 is the additional information provided by the projected angle between the spin and the kick.

In particular, we fully resolve the intrinsic orientation of the system. For 3C 186 to be compatible with a recoiling BH, we find that the spin, the kick, and the line of sight must be close to collinear. The BH spin and kick thus appear perpendicular to each other only because of a prominent projection effect, while in reality they are actually almost parallel. Assuming broad uninformative priors, we show that, if 3C 186 is indeed a recoiling BH, it must be a very unusual one (though its identification might be favored by selection biases). In light of our findings, we formulated additional observational tests that could unveil the true nature of 3C 186, ultimately confirming or ruling out the recoiling BH hypothesis. In particular, the spatial and brightness asymmetry of the approaching and receding radio jet can be used to measure the viewing angle of the jet, which is presumably being launched along the direction of the BH spin.

Our paper is organized as follows. In Sect. 2 we review the current observations of 3C 186. In Sect. 3 we present our statistical modeling and the resulting constraints on the properties of 3C 186. In Sect. 4 we discuss the astrophysical implications of our results and present potential observational tests using radio data. In Sect. 5 we present some concluding remarks. Some additional sanity checks are provided in Appendix A.

2. Quasar 3C 186

The radio-loud quasar 3C 186 is located in a well-studied cluster of galaxies at redshift z = 1.06 (Hewett & Wild 2010). 3C 186 is powered by a BH of mass ∼3 − 6 × 10⁹ M_⊙ inferred using both spectroscopic and photometric data (Chiaberge et al. 2017). Images from the Hubble Space Telescope with the Wide Field Camera 3 and the Advanced Camera for Surveys (Chiaberge et al. 2017; Morishita et al. 2022) show that the AGN is located at a projected separation of 11.1 ± 0.1 kpc away from the photometric center of the host galaxy, coinciding with the kinematic center of its molecular gas (Castignani et al. 2022). Spectroscopic data were obtained with the Hubble Space Telescope Faint Object Spectrograph (Chiaberge et al. 2017), the Sloan Digital Sky Survey (Chiaberge et al. 2017), the Keck/OSIRIS Integral Field Units (Chiaberge et al. 2018) and the Northern Extended Millimeter Array (Castignani et al. 2022). In particular, Chiaberge et al. (2017) reported a velocity shift between the broad and narrow emission lines of −2140 ± 390 km s⁻¹, where the minus sign indicates that the quasar is receding slower than the galaxy.

The presence of both astrometric and spectroscopic offsets drew attention on 3C 186 as a prominent recoiling BH candidate. As discussed at length by Chiaberge et al. (2017), there are other plausible scenarios. AGN disk features such as extreme emitters or winds can explain some of the shifted lines. Regarding the apparent offset, the observed system could be a superposition of two galaxies at different redshifts, where the one hosting the AGN is very compact. While these properties could explain the observations of 3C 186, they also require ad hoc assumptions and/or the presence of additional undetected structure (see Chiaberge et al. 2017 for a thorough discussion).

The radio jet of 3C 186 was observed by the Karl Jansky Very Large Array (VLA) at 8.4 GHz. Data are publicly available at the National Radio Astronomy Observatory archive with project ID AA129 (PI: C.E. Akujor). Earlier, higher-resolution observations from the European Very Long Baseline Interferometry Network (EVN) and the Multi-Element Radio Linked Interferometer Network (MERLIN) at 1.7 GHz have been reported by Spencer et al. (1991). Unfortunately, the original data are not available anymore (Spencer, priv. comm.); subsequent MERLIN observations at a higher frequency, with a comparable resolution, are present in the MERLIN archive and were kindly shared with us (Williams, priv. comm.). Since these observations are less sensitive, and hence produce shallower constraints on the quantities of interest, in the following we make use of the VLA data (which we reduced; see Sect. 4.2) and some summary information from the EVN+MERLIN observations as reported in Spencer et al. (1991).

3. Relativistic modeling

We test whether the observed properties of 3C 186 can be explained by a GW kick as predicted by general relativity. In other words: is there a BH binary that produces a post-merger BH compatible with the observations? A schematic representation of the geometry of the system is provided in Fig. 1.

Fig. 1. Schematic representation of the targeted geometry. The black dot in the center indicates the post-merger BH remnant. The spin (χ) is indicated in orange, the kick velocity (v) is indicated in green, and the line of sight (n) is indicated in blue. Subscripts ∥ and ⊥ indicate projections with respect to the line of sight.

3.1. Bayesian sampling

The properties of the putative pre-merger binary and the direction of the line of sight need to satisfy the following two constraints: First, the kick velocity projected along the line of sight is

$$\begin{array}{c}\hfill v_{\Vert }=-2140\pm 390\mathrm{km}{\mathrm{s}}^{-1},\end{array}$$(1)

as inferred from spectroscopic data (Chiaberge et al. 2017). Second, the angle between the kick and the spin projected onto the plane of the sky is

$$\begin{array}{c}\hfill {\theta }_{v\chi ,\perp }=arccos({\widehat{\mathit{\chi }}}_{\perp }·{\widehat{\mathit{v}}}_{\perp })={90}^{°}\pm {10}^{°},\end{array}$$(2)

as first pointed out by Morishita et al. (2022).

While Morishita et al. (2022) do not explicitly provide an error budget for θ_vχ, ⊥, they report significant digits up to 1°. We thus opted for a conservative uncertainty of 10°. In Appendix A we show that assuming an error on θ_vχ, ⊥ of 1° instead of 10° returns probability distributions that are largely indistinguishable.

We present a Bayesian analysis using a Monte Carlo Markov chain (MCMC) sampler to identify the regions of the parameter space that are compatible with the observations. We assumed a binary progenitor on a quasi-circular orbit characterized by a mass ratio q = m₂/m₁, where m₁ (m₂) is the mass of the heavier (lighter) BH, and dimensionless spin vectors χ_1, 2. They are described by magnitudes χ_1, 2 ∈ [0, 1), polar angles θ_1, 2 between the spins and the orbital angular momentum, and azimuthal angles ϕ_1, 2 measured in the orbital plane. In the following, we often report the difference Δϕ = ϕ₁ − ϕ₂ because it tracks the occurrence of large recoils (Brügmann et al. 2008). We predicted the spin and kick of the post-merger BH using the machine-learning model NRSur7dq4Remnant (Varma et al. 2019a,b), which is trained on numerical-relativity simulations of merging BHs. We adopted their conventions and defined the spin angles at a time t = −100 G(m₁ + m₂)/c³ before BH merger.

We sampled the parameter space using EMCEE (Foreman-Mackey et al. 2013). We used a bivariate Gaussian likelihood,

$$\begin{array}{c}\hfill \mathcal{L}(d|\mathit{x})=\mathcal{N}(\mathit{x}|\mathit{\mu },\mathbf{\Sigma }),\end{array}$$(3)

where d stands for the data and

$$\begin{array}{c}\hfill \mathit{x}=\left(\begin{array}{c}{v}_{\Vert }\\ {\theta }_{v\chi ,\perp }\end{array}\right).\end{array}$$(4)

The likelihood mean and covariance are

$$\begin{array}{c}\hfill \mathit{\mu }=\left(\begin{array}{c}-2140\mathrm{km}{\mathrm{s}}^{-1}\\ {90}^{°}\end{array}\right),\mathbf{\Sigma }={\left(\begin{array}{cc}390\mathrm{km}{\mathrm{s}}^{-1}& 0\\ 0& {10}^{°}\end{array}\right)}^2.\end{array}$$(5)

We set uninformative priors on all our target parameters. Specifically, we sampled the mass ratio q ∈ [1/6, 1] uniformly (the lower bound was set by the validity of the NRSur7dq4Remnant model and we have verified this does not affect our results), the dimensionless spin magnitudes χ_1, 2 ∈ [0, 1) uniformly, and the spin directions isotropically. The direction of the line of sight is parametrized by two angles θ_n and ϕ_n, and is also assumed to be isotropic. The fact that these angles are measured in the same reference frame of the spins (i.e., θ_n = 0°) implies that the line of sight is parallel to the pre-merger orbital angular momentum.

The MCMC sampler converges to a region of the parameters space of the progenitor binaries that is compatible with 3C 186 being a recoiling remnant. That said, our MCMC run shows a low acceptance rate of 9.6% and a considerable autocorrelation length ∼2 × 10⁴ (where we quote the maximum value across 9 parameters) out of 5 × 10⁵ steps for each walker. This indicates an intrinsic difficulty at locating the maximum of the posterior under our uninformative prior.

Figure 2 shows the posterior distribution of all the sampled parameters. Remarkably, this is bimodal, with two configurations that are compatible with the observed properties. The distinguishing factor is θ_n, the angle between the line of sight and the z-axis that traces the direction of the orbital angular momentum before merger. Medians and 90% credible intervals of some parameters are listed in Table 1, where we separate the two modes of the distribution.

Fig. 2. Joint posterior distributions of the parameters describing the binary progenitors of 3C 186. We show one- and two-dimensional marginal posteriors for the mass ratio (q), spin magnitudes (χ1, 2), spin inclination angles (θ1, 2), the in-plane azimuthal angle of the heavier BH spin (ϕ1), the difference between the azimuthal angles (Δϕ), polar observation angles (θn), and the azimuthal observation angle (ϕn). Contours indicate 90% and 50% credible regions.

We also performed an additional run where we omitted the constraint on θ_vχ, ⊥ and only considered that on v_∥; the results are reported in Appendix A. The posterior is visibly different, and in particular prefers binaries with larger mass ratios; the direction of the line of sight is also affected. This suggests that the jet direction is indeed important to constrain the properties of the system.

3.2. Progenitors of 3C 186

BH remnants recoil because of anisotropic emission of GWs, which in turn is due to asymmetries between the masses and spins of the two merging BHs (Gerosa et al. 2018). Our reconstruction shows evident constraints in the parameter space of the putative BH binary that generated 3C 186. Achieving recoils of 𝒪(1000) km s⁻¹ requires configurations with mass ratios close to unity (q ∼ 1), spin magnitudes close to their maximum values (χ_1, 2 ∼ 1), and spin directions that are sufficiently misaligned with respect to the orbital angular momentum (θ_1, 2 ≳ 50°) and antiparallel to each other (Δϕ ∼ 180°; González et al. 2007; Campanelli et al. 2007; Lousto & Zlochower 2011). Moreover, when the recoil velocity is large, its direction is largely parallel to the pre-merger orbital angular momentum (Gerosa et al. 2018), which in turn tracks the direction of the remnant spin. That is, we expect θ_vχ ∼ 0° or 180°.

The results of our MCMC run imply a kick of ∼1700 km s⁻¹, which is large but relatively far from the maximum kick achievable (∼5000 km s⁻¹; Lousto & Zlochower 2011). For this reason, the sampler favors binaries with somewhat asymmetric masses (q ∼ 0.6) and with a high-spinning primary BH (χ₁ ∼ 0.9). Constraints on the secondary spin magnitude are naturally weaker for asymmetric binaries because the ratio between the spin angular momenta is q²χ₂/χ₁. The spins χ₁ and χ₂ do not lie onto the orbital plane (θ₁ ∼ θ₂ ∼ 80°), which indicates an interplay between the so-called superkicks (which occur for θ_i ∼ 90°; González et al. 2007; Campanelli et al. 2007) and hang-up kicks (which occur for θ_i ∼ 50°; Lousto & Zlochower 2011).

We find a strong correlation between the modes of θ_n and those of ϕ₁ (see Fig. 2). This behavior is set by the strong dependence of BH kicks on the orbital phase at merger. A detailed investigation on this point was first presented by Brügmann et al. (2008) and repeated with updated models by Gerosa et al. (2018). In brief, for systems that can excite large kicks, the center of mass of the binary oscillates during the inspiral because of the frame dragging of the two BHs. The merger halts this oscillation, which implies the resulting kick follows a roughly sinusoidal behavior with the orbital phase at merger, which is in turn degenerate with the azimuthal spin projections. In this case, the two ϕ₁ modes correspond to recoils that are qualitatively similar but in opposite directions (the “top” and the “bottom” of the center-of-mass oscillation), so this angle is highly correlated with θ_n.

3.3. Geometry of 3C 186

The spin and the kick of 3C 186 are either almost parallel (θ_vχ ∼ 0° for the mode with 0° < θ_n < 90°) or almost antiparallel (θ_vχ ∼ 180° for the mode with 90° < θ_n < 180°) to each other (see Table 1). This does not contradict the observation that the angle θ_vχ, ⊥ between the direction of the astrometric offset and that of the radio jet needs to be ∼90°, which is indeed is recovered correctly in our sampling run. We are here breaking some of the intrinsic degeneracies of the system. Our key finding is that object 3C 186 is oriented such that the observer is nearly (anti)aligned with both the remnant spin and kick vectors. The kick (v), the spin (χ), and the line of sight (n) are almost collinear; the projections of the kick and spin vectors into the plane of the sky are some small vectorial components, oriented in a way that they appear perpendicular to each other. As a consequence, we find v ∼ |v_∥|≫v_⊥ (see Table 1).

3.4. A rare occurrence?

As shown in Table 1, the inferred value of the recoil velocity along the line of sight v_∥ ∼ 1700 km s⁻¹ is systematically lower than the injected median of 2140 km s⁻¹ from Eq. (5). The low acceptance rate and large autocorrelation time of the MCMC chains points to an intrinsic difficulty of accommodating the observation of 3C 186 as a recoiling BH, at least under our uninformative prior assumptions. As shown in Fig. 3, the likelihood lies in a region of the parameter space where the weight of the prior is vanishingly small – the two barely overlap. From our prior distribution, the probability of obtaining |v_∥ − 2140 km s⁻¹|< 390 km s⁻¹ is as small as 0.15% while the probability of obtaining |θ_vχ, ⊥ − 90° | < 10° is 5.5%. When considering both constraints, we find a probability of 1.0 × 10⁻³%. The qualitative conclusion one can extract from these numbers is that, if 3C 186 is a recoiling BH, it must be a rare one.

Fig. 3. Prior and posterior distributions for the projection of the kick velocity on the line of sight (v∥) and the angle between the spin and kick projections onto the plane of the sky (θvχ, ⊥). The green distribution shows the uninformative prior described in Sect. 3. The posterior distribution sampled with our MCMC run is shown in red. Darker (lighter) colors correspond to high (low) density regions. The blue scatter point with error bars indicates the observational constraints where our bi-variate Gaussian likelihood is centered.

Selection effects are a crucial caveat to the above statement. Recoiling candidates like 3C 186 are identified in large spectroscopic surveys because of their highly shifted lines. So the very fact that 3C 186 was given attention implies that the putative kick is large, and this itself is a rare occurrence. The typical line-shift threshold adopted is of ∼1000 km s⁻¹ (Eracleous et al. 2012). Selecting only the prior volume with v_∥ > 1000 km s⁻¹, the fraction of systems compatible with the two observational constraints as above increases to ∼0.033%.

More formally, one can quantify the impact of selection biases by computing the Bayes factor between two models in which one does (or not) account for such selection threshold. In particular, we find

$$\begin{array}{c}\hfill \mathcal{B}=\frac{p(3\mathrm{C}186|v_{\Vert }\ge 1000\mathrm{km}{\mathrm{s}}^{-1})}{p(3\mathrm{C}186)}=31.7.\end{array}$$(6)

This quantifies that, indeed, the probability of detecting 3C 186-like systems is enhanced by selection effects in current searches. We note that the computation of ℬ from Eq. (6) can be carried out using an appropriate generalization of the so-called Savage-Dickey density ratio and does not require an additional MCMC run.

4. Astrophysical implications

The modeling discussed above sets additional constrains on the nature of 3C 186 and can be used to identify observational tests of the recoiling hypothesis.

4.1. Merger time

The reconstruction of the geometry of the putative recoiling BH allows us to place a novel constraint on the elapsed time since the merger of the progenitors. The magnitude of the kick velocity projected onto the plane perpendicular to the line of sight is v_⊥ ∼ 200 km s⁻¹ (see Table 1). Astrometric observations measure a displacement of d_⊥ = 11.1 kpc ± 0.1 kpc between the AGN and the isophotal center of the galaxy. We provide a Bayesian measurement of the elapsed time t_merger = d_⊥/v_⊥ by combining our posterior on v_⊥ with samples from a normal distribution p(d_⊥) = 𝒩(d_⊥ | μ = 11.1 kpc, σ = 0.1 kpc), thus assuming a constant velocity after the coalescence. Quoting median and 90% credible interval, we find that the BH merger took place $50.5_{-26.0}^{+141.1}$ Myr ago; the resulting posterior distribution is shown in Fig. 4.

Fig. 4. Elapsed time since the BH merger, as inferred by combining the observed astrometric displacement and the posterior distribution of the projected velocity (v⊥). The red line at ∼40 Myr shows the upper limit of the age of the jet inferred from its synchrotron emission.

We next compared the age of the BH remnant to the age of the radio jet. This can be estimated by modeling the spectrum of the emission from the lobes, under the assumption that it is produced by synchrotron emission from electrons accelerated at the shock where the jet collides with the lobe. The age can then be obtained from the position of the “cooling break” in the spectrum, which corresponds to the peak of the synchrotron emission from electrons whose cooling time is equal to the age of the source. Assuming the (effectively isotropic, turbulence-driven) magnetic field in the lobe to hold a fraction of order unity of the internal energy in the shocked plasma (“equipartition”) Murgia et al. (1999) found t_age, eq ∼ 0.1 Myr. They note, however, that their equipartition assumption is motivated by uncertainties on the magnetic field amplification mechanisms, rather than a physical expectation, and that the estimated age of the source scales with the magnetic field strength, B, as

$$\begin{array}{c}\hfill t_{\mathrm{age}}(B)=t_{\mathrm{age},\mathrm{eq}}{\left(\frac{B}{B_{\mathrm{eq}}}\right)}^{-3/2},\end{array}$$(7)

where B_eq = 0.7 mG (Murgia et al. 1999). More recent advancements in particle-in-cell simulations of magnetized relativistic collisionless shocks indicate that the downstream magnetic field decays relatively rapidly, reaching energy densities that are orders of magnitude below those of equipartition in the region where most of the synchrotron emission is produced (e.g., Sironi et al. 2015). This is also supported by observations of synchrotron emission from relativistic shocks caused by gamma-ray-burst-producing relativistic jets (Santana et al. 2014). On the other hand, if the magnetic field were extremely weak, electron cooling in the lobes would have been dominated by inverse Compton scattering with photons of the cosmic microwave background. This effect is equivalent (Murgia et al. 1999) to that of a magnetic field strength B_CMB(z) = 3.25(1 + z)² μG, which is equal to 13.4 μG at the redshift z = 1.06 of 3C 186. This poses an upper limit t_age ≲ t_age, eq(B_CMB/B_eq)^−3/2 = 3.98 × 10⁷ yr ∼ 40 Myr. This upper limit accommodates values that are compatible with the merger time t_merger derived above (see Fig. 4), thus allowing for an interpretation where the jet switched on soon after the BH merger.

If one instead assumes equipartition (i.e., t_age ∼ 0.1 Myr), our estimates of t_merger implies that 3C 186 is a relatively young AGN that switched on tens of megayears after merger. In the recoiling scenario, such a delay between the coalescence of the BHs and the triggering of accretion has been predicted by, for example, Milosavljević & Phinney (2005). The time-dependent gravitational potential of the binary carves a gap in the circumbinary accretion disk. Past coalescence, angular-momentum transport processes (modeled with an effective viscosity) refill the gap and restart accretion on the local viscous timescale. Milosavljević & Phinney (2005) and Dotti et al. (2006) reported that such a delay in the observer frame should be ∼1.5 Myr, though it is important to stress that those estimates were calibrated on significantly lighter BHs of ∼10⁶ M_⊙. While possible, this explanation requires some degree of fine-tuning or at least additional modeling to explain what recent event launched the jet.

Selection effects could also be relevant here: if accretion started earlier, it would either be fainter (hence harder to observe, possibly not producing strong broad lines) or exhausted by the time the BH reaches the observed position. Assuming the accretion disk dragged by the recoiling BH is (i) well described by a Shakura & Sunyaev (1973) disk and (ii) limited in its outer radius by self-gravity driven fragmentation, then the total mass available for accretion is M_disk ≈ M × H/R, where M is the mass of the BH and H/R ∼ 10⁻³ is the aspect ratio for a typical AGN accretion disk (King & Pringle 2007). A BH accreting close to its Eddington limit would then consume all its fuel on a timescale comparable with the age of the jet, partially motivating its youth.

4.2. Jet orientation

We assumed that the jet axis is aligned with the BH spin, so the jet viewing angle (i.e., the angle between the line of sight and the jet axis) is

$$\begin{array}{c}\hfill {\theta }_{\mathrm{jet}}=\{\begin{array}{cc}{\theta }_{n\chi }\hfill & \mathrm{if}{0}^{°}\le {\theta }_{n\chi }<{90}^{°},\hfill \\ {180}^{°}-{\theta }_{n\chi }\hfill & \mathrm{if}{90}^{°}\le {\theta }_{n\chi }<{180}^{°}.\hfill \end{array}\end{array}$$(8)

From our MCMC we find ${\theta }_{\text{jet}}={7.2}^{°}{}_{-{5.2}^{°}}^{+{9.8}^{°}}$. The jet viewing angle can also be measured using radio data and, in the following, we present two different strategies to do so.

We reduced and analyzed the VLA data of 3C 186 using version 6.5.4 of the CASA software (Bean et al. 2022). In particular, we calibrated the flux density scale using the close-by 0812+367 calibrator as reference and we reconstructed the image using natural weighting and a 50 mas cell size. Constant surface brightness contours of the resulting image are shown in blue in Fig. 5. Blue stars mark the centers of four elliptical Gaussian components that we fitted to the image, indicated by A, B, C, and K. In Fig. 5 we superimpose in black the isocontours of surface brightness obtained from the more sensitive, higher-resolution observation performed with EVN+MERLIN at 1.7 GHz (Spencer et al. 1991). Coordinates in the figure refer to the VLA observation and the EVN+MERLIN contours have been aligned using the position of component A.

Fig. 5. Surface brightness contours of the reconstructed image from archival VLA and EVN+MERLIN observations of 3C 186 at 8.4 GHz and 1.7 GHz, respectively. Blue curves indicate the VLA observation; we show six contours of constant surface brightness at 8.4 GHz, logarithmically spaced between 10% and 80% of the peak surface brightness (33.3 mJy beam−1). Black curves in the background refer to the higher resolution EVN+MERLIN image (Spencer et al. 1991). Blue stars mark the centers of four elliptical Gaussian components, which we label A (lobe of the receding jet), B (lobe of the approaching jet), C (core), and K (knot).

Consistently with Spencer et al. (1991) and given the spectral information reported there, we identify the following features.

• We interpret component C as the jet “core”, whose projected position should be consistent with that of the supermassive BH to within a few milliarcseconds (Sokolovsky et al. 2011).

• The component labeled with K is the “knot” of the jet, marking dissipation of energy by means of, for example, internal shocks (Rees 1978).

• Components A and B are the “lobes”, where the kinetic energy of two oppositely oriented relativistic jets is converted into internal energy due to a strong shock caused by collision of the jet particle with the interstellar medium (Blandford & Rees 1974).

Given the absence of visible emission between the core and lobe A, we interpret lobe A to be associated with the receding jet (or “counterjet”) and lobe B to be associated with the approaching jet.

4.2.1. Approaching and receding jet brightness ratio

The first method for estimating θ_jet is based on the observed brightness ratio ℛ between the approaching and receding jets. Assuming Doppler boosting of otherwise identical jets, one has (Blandford & Königl 1979; Scheuer & Readhead 1979)

$$\begin{array}{c}\hfill \mathcal{R}=\frac{F_{\nu ,\mathrm{a}}}{F_{\nu ,\mathrm{r}}}={\left(\frac{1+\beta cos{\theta }_{\mathrm{jet}}}{1-\beta cos{\theta }_{\mathrm{jet}}}\right)}^{2+\alpha },\end{array}$$(9)

where F_ν, a indicates the flux density of the approaching jet, which we identify with that of component K, and F_ν, r is the flux density of the receding jet, which is undetected. The spectral index α is defined such that F_ν ∝ ν^−α for both jets; for component K Spencer et al. (1991) reports α = 1.2. The jet bulk velocity in units of the speed of light is given by β = (1 − γ⁻²)^1/2, where γ is the associated Lorentz factor.

The posterior probability on the viewing angle and the Lorentz factor given the surface brightness data, d, from the image is given by

$$\begin{array}{c}\hfill p({\theta }_{\mathrm{jet}},\gamma |\mathit{d})=\frac{p(\mathit{d}|{\theta }_{\mathrm{jet}},\gamma )\pi ({\theta }_{\mathrm{jet}})\pi (\gamma )}{p(\mathit{d})},\end{array}$$(10)

where π(θ_jet) and π(γ) are prior probabilities and the likelihood can be written as

$$\begin{array}{cc}\hfill p(\mathit{d}|{\theta }_{\mathrm{jet}},\gamma )& ={\displaystyle \int \int p}(\mathit{d}|F_{\nu ,\mathrm{a}},F_{\nu ,\mathrm{r}})p(F_{\nu ,\mathrm{a}},F_{\nu ,\mathrm{r}}|{\theta }_{\mathrm{jet}},\gamma )\mathrm{d}{F}_{\nu ,\mathrm{r}}\mathrm{d}{F}_{\nu ,\mathrm{a}}\hfill \\ \hfill & ={\displaystyle \int \int p}(\mathit{d}|F_{\nu ,\mathrm{a}},F_{\nu ,\mathrm{r}})\delta [F_{\nu ,\mathrm{a}}\hfill \\ \hfill & -\mathcal{R}({\theta }_{\mathrm{jet}},\gamma )F_{\nu ,\mathrm{r}}]\pi (F_{\nu ,\mathrm{r}})\mathrm{d}{F}_{\nu ,\mathrm{r}}\mathrm{d}{F}_{\nu ,\mathrm{a}},\hfill \end{array}$$(11)

where the second equality follows from θ_jet and γ being uninformative with respect to either of the flux densities alone, while informing their ratio through Eq. (9). The likelihood p(d | F_ν, a, F_ν, r) can be expressed as

$$\begin{array}{c}\hfill p(\mathit{d}|F_{\nu ,\mathrm{a}},F_{\nu ,\mathrm{r}})=\frac{p({\mathit{F}}_{\nu ,\mathrm{r}}|\mathit{d})p({\mathit{F}}_{\nu ,\mathrm{a}}|\mathit{d})}{\pi (F_{\nu ,\mathrm{a}})\pi (F_{\nu ,\mathrm{r}})}p(\mathit{d}).\end{array}$$(12)

Combining Eqs. (10)–(12) and carrying out the integral on F_ν, a yields

$$\begin{array}{c}\hfill p({\theta }_{\mathrm{jet}},\gamma |\mathit{d})=\pi ({\theta }_{\mathrm{jet}})\pi (\gamma ){\displaystyle \int {\mathcal{L}}_{F_{\nu ,\mathrm{a}}}}\left[\mathcal{R}({\theta }_{\mathrm{jet}},\gamma )F_{\nu ,\mathrm{r}}\right]p(F_{\nu ,\mathrm{r}}|\mathit{d})\mathrm{d}{F}_{\nu ,\mathrm{r}},\end{array}$$(13)

where we have defined

$$\begin{array}{c}\hfill {\mathcal{L}}_{F_{\nu ,\mathrm{a}}}(x)={\frac{p(F_{\nu ,\mathrm{a}}|\mathit{d})}{\pi (F_{\nu ,\mathrm{a}})}|}_{F_{\nu ,\mathrm{a}}=x}.\end{array}$$(14)

In words, Eq. (13) is the joint posterior probability of the viewing angle and the bulk Lorentz factor informed by the approaching and receding jet flux density measurements. We can then marginalize this expression over γ to obtain the posterior probability p(θ_jet|d). In practice, we approximated these integrals using Monte Carlo sums.

Our Gaussian fit of the VLA radio data yields a total flux density for the K component of F_ν, K = 21.4 ± 1.3 mJy, where we report one-sigma confidence intervals. Based on this, we approximate

$$\begin{array}{c}\hfill {\mathcal{L}}_{F_{\nu ,\mathrm{a}}}(x)\propto \mathcal{N}(x|\mu =21.4\mathrm{mJy},\sigma =1.3\mathrm{mJy}).\end{array}$$(15)

At the putative position of the receding jet, that is, at the mirror position of knot K with respect to the core C, the surface brightness in our VLA image is I_ν, r ∼ 0.5 mJy beam⁻¹. Assuming the image noise to be reasonably well modeled as a zero-mean Gaussian with standard deviation equal to the measured root-mean-square noise σ_rms = 0.6 mJy beam⁻¹ in the image, assuming the receding jet to be un-resolved, and adopting a uniform prior on its flux density, one gets

$$\begin{array}{c}\hfill p(F_{\nu ,\mathrm{r}}|\mathit{d})\sim {\left(2\pi {\sigma }_{\mathrm{rms}}^2\right)}^{-1/2}exp[-\frac{1}{2}{\left(\frac{F_{\nu ,\mathrm{r}}/{\theta }_{\mathrm{b}}^2-I_{\nu ,\mathrm{r}}}{{\sigma }_{\mathrm{rms}}}\right)}^2],\end{array}$$(16)

where ${\theta }_{\text{b}}^2$ is the beam size, which does not need be specified as it simplifies. We adopted a mildly informative, scale-free prior probability on the bulk jet Lorentz factor π(γ) that is uniform in log between 2 and 30. This range accommodates values of the Lorentz factor estimated from superluminal motion in AGN jets (Vermeulen & Cohen 1994).

The red distribution in Fig. 6 shows the resulting posterior on the jet viewing angle θ_jet. Since the measurement of the receding jet flux density is not constraining, the result favors relatively large viewing angles owing to the fact that these are sufficient to produce mild flux-density contrasts. We find that ∼42% of the posterior weight lies in the region where θ_jet > 60°, which is actually disfavored since the broad line region is unobscured (e.g., Baker & Hunstead 1995). For comparison, the blue distribution shows an analogous result obtained the higher-resolution observations by Spencer et al. (1991). In particular, we estimated F_ν, K = 80 ± 0.15 mJy, where the central value is reported in the cited paper and the error is set equal to the image noise rms, and F_ν, K = 0.15 ± 0.15 mJy for the receding jet, as the corresponding feature is not visible in the image. While higher-sensitivity data result in a tighter constraint on the viewing angle, neither of the two datasets allows for a precise determination.

Fig. 6. Posterior distribution of the jet viewing angle (θjet). The red (blue) curve shows results from the approaching and receding jet brightness method (Sect. 4.2.1) using data from archival VLA (EVN+MERLIN) observations. The yellow histogram shows results from the lobe asymmetry method (Sect. 4.2.2), adopting a log-uniform prior on the source age in the range 0.1 ≤ tage/Myr ≤ 40. The green histogram shows the posterior on θjet from our relativistic model (Sect. 3). The gray shaded region corresponds to viewing angles that are disfavored by the observations of broad lines (BLR) in the optical spectrum of the core.

4.2.2. Lobe asymmetry

Our second method for estimating θ_jet relies on the fact that the receding lobe (component A) is located at a slightly larger distance from Earth than the approaching lobe (component B). The gap between the two must have been increasing over time as the lobes were pushed farther out by the jets. Photons reaching us from A have thus been traveling for longer than those from B, and this produces an asymmetry in the apparent lobe separations from the core, which is in principle measurable. We assumed the two lobes moved with the same average speed, β_L (in units of c), from their initial position, which was presumably close to the core. We also assumed that the two jets have similar properties and that the background inter-galactic medium is the same on both sides of the core.

The ratio of the lobe separations s_A and s_B from the core is given by Scheuer (1995)

$$\begin{array}{c}\hfill \frac{s_{\mathrm{B}}}{s_{\mathrm{A}}}=\frac{1+{\beta }_{\mathrm{L}}cos{\theta }_{\mathrm{jet}}}{1-{\beta }_{\mathrm{L}}cos{\theta }_{\mathrm{jet}}},\end{array}$$(17)

so

$$\begin{array}{c}\hfill {\beta }_{\mathrm{L}}cos{\theta }_{\mathrm{jet}}=\frac{s_{\mathrm{B}}/s_{\mathrm{A}}-1}{s_{\mathrm{B}}/s_{\mathrm{A}}+1}·\end{array}$$(18)

The apparent speed of the lobes is

$$\begin{array}{c}\hfill {\beta }_{\mathrm{app},\mathrm{A}}=\frac{{\beta }_{\mathrm{L}}sin{\theta }_{\mathrm{jet}}}{1+{\beta }_{\mathrm{L}}cos{\theta }_{\mathrm{jet}}},{\beta }_{\mathrm{app},\mathrm{B}}=\frac{{\beta }_{\mathrm{L}}sin{\theta }_{\mathrm{jet}}}{1-{\beta }_{\mathrm{L}}cos{\theta }_{\mathrm{jet}}},\end{array}$$(19)

such that

$$\begin{array}{c}\hfill {\beta }_{\mathrm{app},\mathrm{rel}}=2tan{\theta }_{\mathrm{jet}}\frac{{\beta }_{\mathrm{L}}cos{\theta }_{\mathrm{jet}}}{1-{({\beta }_{\mathrm{L}}cos{\theta }_{\mathrm{jet}})}^2}\end{array}$$(20)

is their apparent relative speed. After a time t_age (as measured in the rest frame of the core), the apparent projected physical separation between the lobes is s_ABd_A = β_app, relct_age where s_AB is the separation between the lobes and d_A = 35.38 kpc arcsec⁻¹ is the angular diameter distance at the redshift of 3C 186, assuming a flat Friedmann-Lemaître-Robertson-Walker universe with cosmological parameters from Planck Collaboration XIII (2016). Putting the pieces together, we find

$$\begin{array}{c}\hfill {\beta }_{\mathrm{L}}cos{\theta }_{\mathrm{jet}}=\frac{c{t}_{\mathrm{age}}tan{\theta }_{\mathrm{jet}}}{s_{\mathrm{AB}}{d}_{\mathrm{A}}}[\sqrt{1+{\left(\frac{s_{\mathrm{AB}}{d}_{\mathrm{A}}}{c{t}_{\mathrm{age}}tan{\theta }_{\mathrm{jet}}}\right)}^2}-1]·\end{array}$$(21)

Combining Eqs. (18) and (21) results in an implicit equation that relates the viewing angle θ_jet to quantities that can be measured from the image, namely s_A, s_B, and s_AB, as well as the source age t_age. For small degrees of asymmetry s_B/s_A → 1, the expression can be simplified and solved analytically, leading to

$$\begin{array}{c}\hfill tan{\theta }_{\mathrm{jet}}\sim \left(\frac{s_{\mathrm{AB}}{d}_{\mathrm{A}}}{2c{t}_{\mathrm{age}}}\right)\frac{s_{\mathrm{B}}/s_{\mathrm{A}}+1}{s_{\mathrm{B}}/s_{\mathrm{A}}-1},\end{array}$$(22)

but the solution to the full equation can be easily obtained numerically. The results reported below make use of the numerical solution.

Using VLA data, we measure s_A = 1.023 ± 0.011 arcsec, s_B = 1.157 ± 0.025 arcsec, and s_AB = 2.178 ± 0.023 arcsec, where we quote means and 1σ intervals. We assumed a prior on t_age that is uniform in log between 0.1 Myr and 40 Myr, which accommodates the uncertainties detailed in Sect. 4.1. We thus estimated the jet viewing angle θ_jet by sampling s_A, s_B, and s_AB from Gaussian distributions (with means and standard deviations corresponding to the estimates reported above), sampling t_age from its prior, and evaluating the numerical solution for θ_jet described above. The resulting distribution is shown in orange in Fig. 6. Because of the small degree of lobe asymmetry observed and the “old ages” of 𝒪(10) Myr allowed by the prior (hence relatively low lobe speeds), this estimate favors small viewing angles (θ_v < 69.9° at 90% confidence), which are in good agreement with the recoiling scenario.

5. Conclusions

Black hole merger kicks are a direct consequence of linear-momentum dissipation in strong-field gravity. 3C 186 – a quasar exhibiting both spatial and velocity displacements from its host galaxy – is arguably the most promising astrophysical candidate to date (Chiaberge et al. 2017). The source also shows a prominent radio jet with several resolved features (Spencer et al. 1991).

When modeled as a recoiling BH in general relativity, we find that there exists a region in the parameter space of the possible progenitor binaries that explains all the pieces of observational evidence at our disposal, namely the measured line shift (Chiaberge et al. 2017) and the projected direction of the jet (Morishita et al. 2022). The previous analysis by Lousto et al. (2017) considered progenitor binaries resulting in a remnant recoiling with an (intrinsic, not projected) velocity v > 2000 km s⁻¹. Using their model with a uniform prior and quoting 1σ errors, they report a preference for progenitor binaries with mass ratios $q=0.{49}_{-0.18}^{+0.26}$ and spins that are close to unity $({\chi }_1=0.{99}_{-0.46}^{+0},{\chi }_2=0.{99}_{-0.26}^{+0})$ and misaligned with respect to the orbital angular momentum. As a consequence of the large recoil, they also infer that the imparted kick is nearly collinear with the orbital angular momentum, which is expected as their study predates the consideration on the jet direction by Morishita et al. (2022). When taking projection effects into account, we find $v={1759}_{-643}^{+581}$ km s⁻¹, which indicates that the threshold imposed by Lousto et al. (2017) might be too restrictive. For a comparison with the above values, our best estimates are $q=0.{57}_{-0.29}^{+0.37}$, ${\chi }_1=0.{89}_{-0.41}^{+0.10}$, and ${\chi }_2=0.{60}_{-0.54}^{+0.37}$ at 90% credibility.

In particular, we fully reconstruct the geometry of the system, showing that both the spin and the kick vectors must be pointing toward or away from us. However, the recoiling BH scenario requires some amount of fine-tuning, at least under uninformative prior assumptions. We conclude that, if 3C 186 is indeed a recoiling BH, it must be a rare one.

Such a low probability of obtaining a configuration similar to that of 3C 186 implies the existence of a much larger population of recoiling BHs with projected velocities v_∥ ≳ 1000 km s⁻¹ and redshifts z ≈ 1 within the Sloan Digital Sky Survey footprints. They should be the remnants of close-to-equal-mass binaries of ≳10⁹ M_⊙ that merged within the last ∼50 Myr. If all these systems were to be accreting, the expected number would be far larger than the tens of systems identified within the Sloan Digital Sky Survey presenting large velocity shifts between the broad and narrow emission lines (Tsalmantza et al. 2011; Eracleous et al. 2012). An important caveat to this statement is that the majority of such high-velocity recoiling BHs could be inactive. For future work, our modeling can be used to compute the expected merger rate of such 3C 186-like objects and check whether it is compatible with current GW constraints from pulsar timing arrays (EPTA Collaboration and InPTA Collaboration 2024; Agazie et al. 2023) as well as the rate of massive galaxy mergers (López-Sanjuan et al. 2015).

We then attempted an observational verification of our findings using archival radio data. First, our relativistic model indicates that the BH merger must have happened 𝒪(10) Myr ago. We believe this is compatible with the age inferred from the emission spectrum of the radio lobe, provided one appropriately rescales the nominal age (Murgia et al. 1999) with the expected magnetic field. Second, assuming the radio jet is launched along the direction of the BH spin, our modeling predicts that the viewing angle of the jet must be small. We tested this using both the spatial and the flux asymmetry in the approaching and receding jets, and found that our constraining power is ultimately limited by the quality of the available archival data. Deeper radio observations of 3C 186 might shed light on its nature as a recoiling supermassive BH.

Acknowledgments

We thank Vishal Baibhav, Emanuele Berti, Marco Chiaberge, David Izquierdo-Villalba, Colin Norman, Alberto Sesana, Ralph Spencer, David Williams and Georg Herzog for discussions. M.B. and D.G. are supported by ERC Starting Grant No. 945155–GWmining, Cariplo Foundation Grant No. 2021-0555, MUR PRIN Grant No. 2022-Z9X4XS, and the ICSC National Research Centre funded by NextGenerationEU. D.G. is supported by MSCA Fellowships No. 101064542–StochRewind and No. 101149270–ProtoBH. O.S.S. acknowledges support by the PRIN-INAF grant 1.05.23.04.04. Computational work was performed at CINECA with allocations through INFN and Bicocca.

References

Appendix A: Additional runs

We present results from some additional MCMC runs, extending what is presented in Sect. 3. In Fig. A.1 we repeat our analysis assuming the standard deviation of θ_vχ, ⊥ is 1° instead of the more conservative value 10° used in the many body of the paper. The results are indistinguishable. In Fig. A.2 we ignore the constraint on the projected spin-kick angle (θ_vχ, ⊥ = 90° ±10°) and only consider that on the projected kick velocity (v_∥ = −2140 ± 390 km/s). In this case the posterior is different, with the most affected parameters being the mass ratio (q), the angle between the in-plane spin projection (Δϕ), and the orientation angles (θ_n, ϕ_n).

Fig. A.1. Joint posterior distributions of the parameters describing the binary progenitors of 3C 186 assuming a 10° (blue) or 1° (red) error on θvχ, ⊥. The blue distribution is the same as in Fig. 2.

Fig. A.2. Joint posterior distributions of the parameters describing the binary progenitors of 3C 186 with (blue) and without (gold) the constraint on θvχ, ⊥. The blue distribution is the same as in Fig. 2.

All Tables

All Figures

	Fig. 1. Schematic representation of the targeted geometry. The black dot in the center indicates the post-merger BH remnant. The spin (χ) is indicated in orange, the kick velocity (v) is indicated in green, and the line of sight (n) is indicated in blue. Subscripts ∥ and ⊥ indicate projections with respect to the line of sight.
In the text

Fig. 2. Joint posterior distributions of the parameters describing the binary progenitors of 3C 186. We show one- and two-dimensional marginal posteriors for the mass ratio (q), spin magnitudes (χ1, 2), spin inclination angles (θ1, 2), the in-plane azimuthal angle of the heavier BH spin (ϕ1), the difference between the azimuthal angles (Δϕ), polar observation angles (θn), and the azimuthal observation angle (ϕn). Contours indicate 90% and 50% credible regions.

In the text

Fig. 3. Prior and posterior distributions for the projection of the kick velocity on the line of sight (v∥) and the angle between the spin and kick projections onto the plane of the sky (θvχ, ⊥). The green distribution shows the uninformative prior described in Sect. 3. The posterior distribution sampled with our MCMC run is shown in red. Darker (lighter) colors correspond to high (low) density regions. The blue scatter point with error bars indicates the observational constraints where our bi-variate Gaussian likelihood is centered.

In the text

	Fig. 4. Elapsed time since the BH merger, as inferred by combining the observed astrometric displacement and the posterior distribution of the projected velocity (v⊥). The red line at ∼40 Myr shows the upper limit of the age of the jet inferred from its synchrotron emission.
In the text

Fig. 5. Surface brightness contours of the reconstructed image from archival VLA and EVN+MERLIN observations of 3C 186 at 8.4 GHz and 1.7 GHz, respectively. Blue curves indicate the VLA observation; we show six contours of constant surface brightness at 8.4 GHz, logarithmically spaced between 10% and 80% of the peak surface brightness (33.3 mJy beam−1). Black curves in the background refer to the higher resolution EVN+MERLIN image (Spencer et al. 1991). Blue stars mark the centers of four elliptical Gaussian components, which we label A (lobe of the receding jet), B (lobe of the approaching jet), C (core), and K (knot).

In the text

Fig. 6. Posterior distribution of the jet viewing angle (θjet). The red (blue) curve shows results from the approaching and receding jet brightness method (Sect. 4.2.1) using data from archival VLA (EVN+MERLIN) observations. The yellow histogram shows results from the lobe asymmetry method (Sect. 4.2.2), adopting a log-uniform prior on the source age in the range 0.1 ≤ tage/Myr ≤ 40. The green histogram shows the posterior on θjet from our relativistic model (Sect. 3). The gray shaded region corresponds to viewing angles that are disfavored by the observations of broad lines (BLR) in the optical spectrum of the core.

In the text

	Fig. A.1. Joint posterior distributions of the parameters describing the binary progenitors of 3C 186 assuming a 10° (blue) or 1° (red) error on θvχ, ⊥. The blue distribution is the same as in Fig. 2.
In the text

	Fig. A.2. Joint posterior distributions of the parameters describing the binary progenitors of 3C 186 with (blue) and without (gold) the constraint on θvχ, ⊥. The blue distribution is the same as in Fig. 2.
In the text