The Host Galaxy of the Recoiling Black Hole Candidate in 3C 186: An Old Major Merger Remnant at the Center of a z = 1 Cluster

3C 186 was observed in a General Observer program (PID 15254; PI. M. Chiaberge) during the HST Cycle 25 with WFC3-IR and Advanced Camera for Surveys (ACS). Two filters, F110W and F160W, are used in the WFC3 observations and one filter, F606W, in the ACS observations (Table 1). The observing position angle was carefully designed during the phase II process so as to minimize possible contamination of the host galaxy by stellar spikes from the quasar and surrounding objects. To avoid saturation of a few pixels in the quasar core, we started with one short exposure (∼30 s) in each of the five ACS observations (one per orbit). Then, the following exposures were taken at the same location and rotation angle of the first dither point in each orbit, to ensure an identical point-spread function (PSF) profile over all exposures.

To reduce imaging data, we follow the procedure presented in Hilbert et al. (2016). For ACS data, we begin with files in the flc format. These images are corrected for charge transfer efficiency (CTE) effects through the CTE correction algorithm (Anderson & Bedin 2010). For WFC3-IR data, we begin with files in the flt format, because the multiple nondestructive detector readouts within each observation allowed for easy cosmic-ray identification and removal within the calwf3 data reduction pipeline. However, pixels subjected to high flux levels in preceding observations are affected by persistence and need to be corrected. To remove any persistence signal present in our data, we begin by retrieving the persistence masks and persistence-corrected flt files of our observations from MAST. These files have had the persistence signal modeled and subtracted from them following the model described in Section 8.3 of version 4.0 of the WFC3 Data Handbook, and are therefore different from the standard flt files available in the archive. Before proceeding with the data reduction, we also manually mask the satellite trail present in one of our images.

Since our focus in this study is the host galaxy of 3C 186, we attempt to subtract the light from the bright quasar nucleus in each distortion-uncorrected image (those in flt/flc formats) before the geometric correction and final stacking steps (see Section 2.2).

The next step in our data reduction is to remove the geometric distortion from all of the PSF-subtracted flc/flt files and to combine images into a final image for each filter. We use Astrodrizzle to accomplish both of these steps.

For each filter, we combine the PSF-subtracted images into a final image with a pixel scale of 0045 pixel⁻¹. We find a value of "5.5 5.0" for the driz_cr_snr parameter results in good cosmic-ray rejection during the stacking process. We also calculate an optimal final_pixfrac value of 0.75 for our data and final pixel scale.

We then use Tweakreg (also part of the Drizzlepac software package) to align this final image to the same world coordinate system present in the corresponding ACS drizzled image. At this point, it is possible to overlay the ACS and WFC3 images for a particular object and compare the morphology and brightness in the two observation bands.

The resulting final image is shown in Figure 1. The PSF-subtracted image reveals a clear spatial offset between the quasar and the host galaxy (Section 3.1). Our deep images also reveal faint outskirts of the host, out to ∼60 kpc.

Zoom In Zoom Out Reset image size

In order to subtract the PSF component without oversubtracting the host light component, we fit the observed light profile with PSF+single Sérsic profile using GALFIT (Peng et al. 2002). In an effort to achieve the best model fits, we experiment with the GALFIT inputs. We create a custom bad pixel file for each exposure, where we remove all nearby astrophysical sources as well as bad pixels and cosmic rays. By supplying this file to GALFIT, we reduce the chances that the final fit is contaminated by nearby sources.

The PSFs used for our PSF subtractions are created using TinyTim software (Krist et al. 2011). We use TinyTim PSFs due to a lack of unsaturated, isolated PSFs with high signal-to-noise ratio (S/N) in our data to use as models. In order to create realistic PSFs, we used characteristics from our observations when creating the model PSFs. We begin by finding the mean secondary mirror focus offset for the date and time of our observation. This information is collected from the HST focus model website. ¹²

Other information used as input to TinyTim includes the HST instrument and filter, the target's pixel location on the detector, the amount of jitter (0 for this work), and the size in arcseconds of the output model PSF (7''). Finally, we also provide a quasar spectrum, in order for TinyTim to determine wavelength-dependent PSF characteristics. For simplicity, we used the composite quasar near-IR spectrum from Glikman et al. (2007) redshifted to z ∼ 1. With this information in hand for each exposure, we run all three stages of TinyTim. The final output is one geometrically distorted model PSF corresponding to each flc/flt image.

We also fit the data using a PSF at the nominal pixel scale, as well as a supersampled PSF. We find that the fit with the supersampled PSF tends to oversubtract the source flux, and therefore keep to the nominal pixel scale for our final fittings. We also fit the data with TinyTim PSFs from a range of focus values, in order to examine any differences. The best modeling and subtraction came from values that match the focus values indicated by the focus model website. Finally, as the data being fit are in units of electrons per second, rather than the electrons that GALFIT assumes, we also experiment with fits where we provide the ERR array of the flt/flc file to GALFIT as a sigma image. We find in this case that the resulting fits again tend to oversubtract the source's flux. With our final set of GALFIT parameters and inputs, we obtain fits with reduced χ² between 0.75 and 1.5. The best-fit PSF component is then subtracted from each of the images before the images are processed for geometric correction and final stacking with Astrodrizzle.

We estimate the center of mass of the host galaxy by using F160W, which is the deepest among our images here. It is noted that the light distribution in the F160W image is considered to trace the stellar component of the host well, because at the redshift of 3C 186 the filter corresponds to rest-frame ∼0.8 μm and is free from strong emission lines (see also Section 3.2).

While the extended light from the host galaxy is clearly visible in the PSF-subtracted image (Figure 1), it is still challenging to accurately estimate its light-weighted center (LWC) due to the fact that fore/background objects are superposed on the host. Especially, since 3C 186 is in a crowded field, the latter remains as a critical issue because flux contribution from other sources can easily affect the estimate.

To eliminate such flux contamination, we first manually mask regions that are likely to belong to surrounding objects (including the compact star-forming blob discussed in Section 3.3) and residuals at the position of the quasar, by placing circular apertures of arbitrary sizes in the PSF-subtracted F160W image. By doing so, we ensure to mask all possible contaminants, whereas automated software (e.g., SExtractor; Bertin & Arnouts 1996) may not necessarily recognize all objects in crowded regions or might mistakenly mask part of the host. The image with those aperture masks is shown in the left panel of Figure 2. To mask the residuals at the location of the quasar, we place an aperture with radius of r = 09.

Zoom In Zoom Out Reset image size

The next step is to fill those masked regions. This process is necessary because leaving those masked pixels could result in a biased estimate of LWC depending on the choice of position and size used of each circular mask. For each pixel of the masked regions, we use a Metropolis–Hastings-like algorithm to randomly select one of the neighboring pixels that are not masked, to fill in the target pixel. This filling process is as follows.

• 1.  
  For the pixel of interest in a masked region, assign a radius (r_rand) randomly taken from the range $[0,{r}_{\max }]$, where we here set ${r}_{\max }=50$ pixels.
• 2.  
  Compare the probability of r_rand, assigned by $p(r)=\exp (-r/{r}_{\max })$, with a random float number, p_r , drawn from the range [0, 1].
• 3.  
  If p(r_rand) > p_r , then select a random pixel from unmasked regions at r_rand < r < r + Δr, where we set Δr = 2 pixels; otherwise return to step 1 and repeat the process until it passes the criterion.
• 4.  
  Repeat steps 1–3 for all masked pixels.
• 5.  
  Repeat steps 1–4 n_iter times, to get n_iter images with the masked regions filled, respectively.

We here set n_iter = 100. By adopting an exponential function for p(r), the filling flux is likely to be chosen from one of the nearby pixels, making the reproduced light profile contiguous, which is reasonable given the primary purpose here. One of the realizations is shown in Figure 2 as an example. While several pixels can occasionally have a value distinct from the average value of surrounding pixels, contributions from such pixels to the determination of LWC are smoothed out after n_iter realizations.

At each realization of n_iter, we calculate the LWC from the following equations:

$$\begin{eqnarray}\begin{array}{rcl}\overline{x} & = & \displaystyle \frac{\sum {{xf}}_{160}(x,y)/{e}_{160}{\left(x,y\right)}^{2}}{\sum {f}_{160}(x,y)/{e}_{160}{\left(x,y\right)}^{2}},\\ \overline{y} & = & \displaystyle \frac{\sum {{yf}}_{160}(x,y)/{e}_{160}{\left(x,y\right)}^{2}}{\sum {f}_{160}(x,y)/{e}_{160}{\left(x,y\right)}^{2}}\end{array}\end{eqnarray} \tag{ 1 }$$

where f₁₆₀(x, y) and e₁₆₀(x, y) represent the flux value and its associated rms value at the position [x, y], and the sum goes over pixels with f₁₆₀/e₁₆₀ > S/N_limit within ${r}_{\mathrm{lim}}\lt 250$ pixels (∼91 kpc) from the target. We here set S/N_limit = 0.7, which is equivalent to the detection threshold (DETECT_THRESH of SExtractor) used for source detection in the F160W image. We empirically find that the value in general provides reasonable detection out to faint outskirts of sources without being contaminated by background pixels.

We calculate the median LWC over n_iter realizations and find $[x,y]=[-{29.42}_{-0.08}^{+0.07},{7.44}_{-0.05}^{+0.06}]$ (in pixels) with respect to the position of 3C 186 in the image of Figure 2, or 7:44:17.586,+37:53:18.165 in sky coordinates. The 16th/84th percentiles are associated as uncertainties. The result reveals a significant offset of 137 ± 001, or 11.1 ± 0.1 kpc, in the projected distance, which is consistent with the estimate with a shallower F140W image by Chiaberge et al. (2017). The angle from the LWC toward the quasar core, measured from the west in the clockwise direction, is ∼+35°, which lies almost perpendicular to the direction of the radio jet (Figure 3). We will discuss this in Section 4.1.

Zoom In Zoom Out Reset image size

The choice of ${r}_{\mathrm{lim}}$ and the limiting S/N has only a minor effect on the estimate as long as S/N_limit is set to >0.7; having lower limiting S/N would significantly increase the estimated error due to noisy pixels across the entire image. On the other hand, the inferred offset becomes even larger when a higher value for S/N_limit is chosen. As in the middle panel of Figure 2, an ellipse fit to the F160W light profile shows a shift of the centroid in the northeast direction when only the central part of 3C 186 is used. This is due to the light concentration in the outskirts. We discuss the light concentration in relation to merger-induced morphological disturbances, or lack thereof, in Section 4.

Despite the depth that our new images reach, we might still be missing components with very low surface brightness (>25 mag arcsec⁻²). However, we argue that such light components only marginally affect the determined photocenter. Our images here already reveal a substantial amount of stellar mass ¹³ (${\rm{log}}({M}_{\ast }/{M}_{\odot })\sim 11.4$; see Section 3.2). This amount is 10× more than what is found in the diffuse outskirts of low-z elliptical galaxies (e.g., Huang et al. 2018).

We aim to investigate the underlying stellar population in the host. The spatial distribution of stellar age enables us to estimate the timeline since the last merger and associated star formation activity.

To get sufficient S/N for reliable estimates out to faint outskirts of the host, we define subregions of the host by using a Voronoi tessellation method (Cappellari & Copin 2003). We set the minimum S/N to 15, so that every tessellated region has S/N greater than this value. As a result, at ∼6'' (∼48 kpc) from the system LWC we reach down to ∼25 mag arcsec⁻² in F160W. We set the boundary of the host by using the segmentation mask generated by SExtractor with the detection threshold of 0.7, which is consistent with the threshold used for the LWC calculation above. Note that in this analysis, as it is not necessary for our purposes, we do not fill the masked regions defined in Section 3.1 and leave those empty. We show the defined mask and radial distribution of tessellated pixels in Figure 4.

Zoom In Zoom Out Reset image size

The boundary for tessellation defined in the F160W image is then applied to both F110W and F606W images, to collect flux from the same subregions. Each of the F110W and F606W images is convolved to the PSF size of the F160W image beforehand. The convolution kernels were generated by providing median-stacked PSFs to the Python software pypher (Boucaud et al. 2016), in the same way as in Morishita (2021). Fluxes of each subregion are summed and compiled into a catalog for a spectral energy distribution (SED) analysis in the following section.

The flux catalog collected from the tessellated images is provided to the SED fitting software FAST (Kriek et al. 2009). Due to the number of broadband filters available here (three), we are limited to a small number of fitting parameters. Since our primary focus here is to estimate the epoch of the last significant star formation (which could potentially be associated with major merger), we simply assume a single burst that started at T₀ with a length of τ. The schematic of this star formation model is shown in Figure 5.

Zoom In Zoom Out Reset image size

We use the Bruzual & Charlot (2003) stellar population model, with the redshift fixed to the one estimated with the Sloan Digital Sky Survey DR6, z = 1.07 (Hewett & Wild 2010), and the Calzetti et al. (2000) dust attenuation model. It is generally challenging to determine metallicity with only broadband photometry, and thus we fix it to solar. This is reasonable given that galaxies in this mass range are already metal-enriched, as seen at similar or earlier epochs (e.g., Thomas et al. 2010; Onodera et al. 2015; Kriek et al. 2019; Morishita et al. 2019).

In Figure 5, we show examples of the SED results in three different regions: the outer region (∼09 from the LWC; Voronoi ID 4), the inner (∼02 from the LWC; Voronoi ID 13), and the blue blob (see Section 3.3). These examples demonstrate that three data points are sufficient to constrain relative differences in age of different regions, by capturing both the Balmer break by V₆₀₆ − J₁₁₀ and the rest-frame optical slope by J₁₁₀ − H₁₆₀. The best-fit parameters for all tessellated regions are summarized in Table 2.

Table 2. Positions of Tessellated Segments and the Best-fit SED Parameters

Voronoi ID	S/N	ra	${\rm{log}}({M}_{\ast }/{M}_{\odot })$	log(T0/yr)	${\rm{log}}(\tau /{\rm{yr}})$	${\rm{log}}({t}_{Q}/{\rm{yr}})$	AV
		(arcsec)					(mag)
0	19.0	0.59	9.2	8.1	7.8	7.8	1.7
1	22.8	0.80	9.6	8.1	7.0	8.1	1.5
2	26.7	1.91	10.0	8.0	7.2	7.9	1.2
3	23.1	0.86	9.6	8.5	7.8	8.4	1.2
4	23.9	4.40	10.7	9.4	9.8	⋯ b	1.3
5	21.0	0.92	9.4	8.8	8.2	8.7	0.0
6	23.1	2.51	9.9	8.3	7.0	8.3	1.1
7	23.4	0.35	9.2	8.0	7.6	7.8	1.3
8	22.3	7.06	10.6	8.3	7.6	8.2	2.1
9	18.7	0.20	9.0	7.9	7.0	7.8	1.0
10	19.3	0.12	9.2	8.3	7.6	8.2	0.9
11	22.5	1.84	10.0	8.1	8.0	7.4	1.7
12	23.4	1.36	9.7	8.8	9.4	⋯ b	1.6
13	15.7	0.18	8.9	8.0	7.0	8.0	0.9
14	20.3	0.55	9.3	8.6	8.2	8.4	0.0
15	18.2	0.39	9.0	7.9	7.2	7.8	1.0
16	16.4	0.36	9.1	8.5	8.0	8.3	0.4
17	22.1	2.60	10.3	8.1	7.6	7.9	2.1
18	17.9	0.49	9.2	8.0	7.0	8.0	1.2
19	21.5	0.81	9.5	7.9	7.0	7.8	1.6
20	26.0	1.52	9.7	8.0	7.8	7.6	1.8
21	20.8	1.21	9.6	7.9	7.4	7.7	1.8
22	26.1	3.87	10.4	8.1	7.6	7.9	1.6
23	22.3	2.51	10.2	8.3	7.8	8.1	1.9
24	22.4	0.86	9.5	8.5	7.0	8.5	0.0
25	23.8	1.52	9.8	8.4	9.6	⋯ b	1.7
26	18.1	2.86	10.1	8.7	10.0	⋯ b	1.8
27	15.4	4.21	10.3	8.2	7.6	8.1	2.2
28	24.9	2.12	10.1	8.3	8.0	8.0	0.6

Notes.

^a Distance from the light-weighted center. ^b Regions with ongoing star formation (i.e., t_Q = 0).

Download table as:  ASCIITypeset image

A caveat is that with the number of data points used, degeneracy between age and dust is still not completely resolved, which makes inference of absolute values challenging. Our SED fitting results indicate small to moderate attenuation (A_V ≲ 2.0 mag). However, in the case of the 3C 186 host we do not expect significant dust attenuation, because no dust lane-like feature is seen in our HST images—such a feature would be clearly seen in the presence of significant dust attenuation (see, e.g., Skipper & Browne 2018; Hogg et al. 2021, for offset AGN candidates with clear dust lanes). In Figure 6, we show the two-dimensional distribution of J₁₁₀ − H₁₆₀ color of the host. While the color does not necessarily represent the exact dust attenuation due to the degeneracy, from the smooth morphology across the host it is unreasonable to expect significant dust structures.

Zoom In Zoom Out Reset image size

In addition, Podigachoski et al. (2015) reported tentative detection (<3σ) or nondetection of 3C 186 in far-IR bands, at 160, 250, 350, and 500 μm by the Herschel Space Observatory, implying no significant amount of dust is present in the system. Fixing dust attenuation to zero would increase the age estimates, to account for red components; therefore, our age estimates, especially those with nonzero A_V , should be taken as lower limits.

In Figure 7, we show the spatial distribution of each of the derived parameters. Our interest here is in T₀, the time at which primary star formation started, and t_Q , which we introduce to refer to length of time since star formation was truncated (≡T₀ − τ). There is no clear radial trend for T₀, with a median value of ${160}_{-40}^{+200}$ Myr, with the associated errors capturing the 16th/84th percentiles. Similarly, no clear spatial correlation nor radial trend is observed for t_Q . The parameter ranges from ∼30 to ∼500 Myr, with a median value of ${t}_{Q}={90}_{-30}^{+110}$ Myr, except in a few subregions with ongoing star formation (i.e., t_Q ≤ 0). The total mass, derived by summing the best-fit estimate over the Voronoi segments, is found to be ${\rm{log}}({M}_{\ast }/{M}_{\odot })=11.4$.

Zoom In Zoom Out Reset image size

As we highlight in the color-composite image of Figure 1, there is a compact star-forming blob (SFB1 hereafter) located at ∼2'' from the quasar in the northeast direction. This blob was seen in the early data and reported in Chiaberge et al. (2017), but its properties remained unclear. Our SED modeling revealed its age, T₀ ∼ 160 Myr, and stellar mass, ∼8 × 10⁹ M_⊙ (Figure 5). We estimate its star formation rate (SFR), 65 ± 20 M_⊙ yr⁻¹, by averaging over the last 100 Myr of the best-fit star formation history, which is consistent with the upper limit of the whole system derived from nondetection in Herschel data, <80 M_⊙ yr⁻¹ (Podigachoski et al. 2015).

The best-fit SFR converts to a star formation rate surface density of ∼1.2 M_⊙ yr⁻¹ kpc⁻². Considering its stellar surface density ($\mathrm{log}{{\rm{\Sigma }}}_{* }\sim 8.1\,{M}_{\odot }\ {\mathrm{kpc}}^{-2}$), this value is ∼0.8 dex higher than the average value of resolved subregions in star-forming galaxies at z ∼ 1 (e.g., Wuyts et al. 2013), indicating intense star formation activity may still be ongoing in the host despite relatively old stellar populations elsewhere. Such intense star formation is actually seen in a post-merger system (e.g., Barrows et al. 2018) but lasts only for a short timescale of ∼100 Myr. We revisit this in Section 4.

Since the template library used in the SED fitting in Section 3.2 does not include emission lines, we check the fidelity by running another SED fitting code, gsf (Morishita et al. 2019), with emission lines included. While an exponentially declining star formation history is used this time, due to the limitation of the code, it derives a consistent value for star formation rate, 70 ± 10 M_⊙ yr⁻¹.

While the IFU observations presented in Chiaberge et al. (2018) covered SFB1, no line was detected, because the observations were originally designed for the brightness of the quasar and exposure was not deep enough to detect any emission lines in SFB1. Recent observations at NOEMA detected extended CO emission near the flux peak position of SFB1 at the same redshift as the optical narrow lines (Castignani et al. 2022).

As we introduced in Section 1, based on its spectroscopic and imaging properties, 3C 186 was proposed as a clear example of a GW recoiling black hole candidate by Chiaberge et al. (2017). With our deep HST observations, we both confirmed and more accurately determined the presence of a spatial offset between the quasar and the host photocenter. Based on our results (Section 3.1), we see that the direction of the putative GW kick lies perpendicular to the radio jet axis. In Figure 3, we show the 8.44 GHz radio continuum emission overlaid on the central region of the host galaxy. The direction of the kick is at the position angle of ∼+33°, while the radio jet axis has a position angle of ∼−57°. As pointed out by, e.g., Lousto et al. (2012), the maximum velocity of a GW kick (>2000 km s⁻¹) may be more likely achieved if the direction of the kick is aligned with the angular momentum of the orbital plane of the merging black hole binary. Since the radio jet axis likely indicates the orientation of the spin of the merged black hole, these results may impose further constraints on the properties of the progenitor binary, in particular with respect to the spin amplitude and relative orientation. A detailed modeling of the binary black hole system and the GW recoil kick lies beyond the scope of this work.

In this subsection, we discuss the main results of this work, which focuses on the observed properties of the host. This allows us to speculate on the current evolutionary stage of this system. Chiaberge et al. (2017) reported a possible shell or tidal feature in the host light profile and suggested its origin as a major merger that occurred about ∼1 Gyr ago from comparison with numerical simulations. Our F160W image presented here revealed even fainter stellar components out to ∼60 kpc from the LWC of the host. While a slight light concentration is seen in the outermost part in the southeast direction, we did not see clear evidence of disturbed morphology or young star clusters in its outskirts, which would be characteristic of ongoing merger (e.g., Mulia et al. 2015), implying that the system is a relatively old, post-merger remnant. Furthermore, Podigachoski et al. (2015) derived high SFR in other 3C galaxies at similar redshifts, most of which are associated with ongoing mergers, whereas an upper limit was estimated for 3C186 from nondetection in their Herschel data (<80 M_⊙ yr⁻¹). This suggests that 3C 186 is older than other 3C systems in terms of recent star formation too.

To further investigate the current evolutionary phase of 3C 186, we compare the observed properties with hydrodynamical simulations. In Figure 8, we show a time sequence of a merging system from the VELA-Sunrise project (Simons et al. 2019). The VELA-Sunrise project provides a set of mock images of simulated galaxies that were originally taken from the VELA simulation suite (Ceverino et al. 2014; Zolotov et al. 2015). The simulation suite consists of a variety of galaxies at different evolutionary stages, including both isolated and merging systems. While the original simulations do not include specific recipes for black hole mergers, the data sets offer detailed morphological information of, e.g., a galaxy during merger sequences, at the resolution of HST. Here we select VELA04, one of the merging systems available in the VELA suite. Specifically, the data set of VELA04 captures the timeline of pre- to post-merger and allows us to investigate the evolutionary stage of 3C 186 by comparing their morphological properties.

Zoom In Zoom Out Reset image size

The example shown in Figure 8 is a sequence that captures the pre-merger to post-merger phases of two galaxies at z ∼ 2. The sequence exhibits several features that resemble the observed properties of the 3C 186 host. The last frame of Figure 8 captures the elongated (but not disturbed) morphology of the host galaxy and compact star-forming regions in the main body of the host. We also see the presence of young, compact blue regions across the entire sequence. Those blue regions do not appear in more than one simulation frame, implying that their lifetime is less than its time resolution, ≲200 Myr. In the last frame of the sequence in Figure 8, we highlight one such compact blue region. Its specific star formation rate is ∼40 Gyr⁻¹ (cf. ∼8 Gyr⁻¹ of the blue blob in 3C 186). The rest of the stellar populations of the galaxy in the frame are relatively old, with a median age of ∼1 Gyr (cf. ≳200 Myr in 3C 186). While we also notice morphological differences between the two systems (e.g., more clumps and less relaxed), those physical properties are broadly consistent with the observed properties of 3C 186 (Section 3.3).

The spatial offset found in this study is encouraging regarding a typical timescale after a recoiling event. By using the derived spatial offset of the quasar from the LWC (Section 3.1) and the velocity offset derived in Chiaberge et al. (2017), we estimate that the recoil event occurred ≳5 Myr ago, though the exact age depends on the assumed kick angle. Theoretically, it is possible for a recoiling clack hole to remain active for such a length of time. By using the formula in Loeb (2007) with assumed radiative efficiency of = 0.1 as in Chiaberge et al. (2017), the lifetime of an accretion disk for the black hole mass of 3C 186 is estimated to be ∼100 Myr. In fact, Murgia et al. (1999) derived the age of the radio source in 3C 186, ∼0.1 Myr, from its synchrotron spectrum, which implies that the radio AGN turned on after the recoiling event.

In addition, Blecha et al. (2011) presented several examples of simulated recoiling black holes and estimated a typical lifetime of ∼10 Myr (see also Blecha et al. 2016). We note, however, that uncertainties in absolute timescale in those simulations are not negligible. Such uncertainties primarily originate in the resolution element, which is typically much larger than the scale of black hole merger (several kiloparsecs), making the comparison with observations challenging. This is still the case with more recent simulations that record the exact time of the black hole merger (e.g., TNG; Hani et al. 2020). Sub-grid models of black hole binary evolution (e.g., Kelley et al. 2017) show a wide range of black hole merger timescales, from ∼1 Gyr to more than a Hubble time, depending on the model assumptions and the specific conditions in the galactic nuclei.

As we mentioned in Section 1, there is an alternative scenario for the observed spatial offset, where the quasar is hosted by another, fainter galaxy ("galaxy B" for convenience) that is superposed on a galaxy at different redshift ("galaxy A" i.e., the one we have analyzed in this study as the host). However, this scenario is unlikely because (1) we did not confirm any stellar light concentration near the position of the quasar in any filters and (2) there is evidence of narrow absorption at the same redshift as the narrow emission lines.

For (1), assuming the scaling relation between M_* and M_BH (e.g., Kormendy & Ho 2013), stellar mass M_* ∼ 10¹² M_⊙ is expected for the black hole mass of 3C 186, M_BH ∼ (3–6) ×10⁹ M_⊙. Such a mass should be noticeable in our F160W image regardless of its light profile. Chiaberge et al. (2017) tested this with various structural parameters for unseen galaxy B by using galfit (Peng et al. 2002) but concluded that one is not likely not to see such a mass even in their F140W image. The recoiling black hole scenario, which assumes galaxy A (${\rm{log}}({M}_{\ast }/{M}_{\odot })\sim 11.4;$ Section 3.2) as the host, offers more reasonable agreement in terms of the scaling relation for the stellar and black hole masses observed. To further investigate this scenario, we repeat the structural analyses by using our deep F160W image. We provide the original image with the quasar component unsubtracted as the input image to galfit. We fit the light profile of the 3C 186 system with one PSF component at the position of the quasar, one Sérsic component fixed at the LWC position derived above, and one additional Sérsic component (i.e., galaxy B) at the quasar position. We fix the position angle and axis ratio of the host to 36° and 0.37 (Chiaberge et al. 2017), respectively, and leave all other parameters free. We find that the second Sérsic component does not improve the fit. The best-fit effective radius of this additional component is unphysically small, 54 mas. In addition, the total flux of this component remains at ∼7% of the host. This converts to a stellar mass of ∼1.8 × 10¹⁰ M_⊙ by adopting the same mass-to-light ratio as the host, which is significantly undermassive for the black hole mass of 3C 186. We therefore conclude that this additional component, if not an artifact such as flux residuals of the PSF component, is unlikely to be the unseen host galaxy of the quasar.

For (2), Chiaberge et al. (2017) found narrow absorption features in both Lyα and [C iv] lines at the same redshift as the narrow emission lines ([O ii] and [Ne iii]). This indicates the presence of gas associated with galaxy A and thus secures its redshift (see also Castignani et al. 2022, who identified CO emission at the same redshift as the narrow lines within the host). Determining the systemic redshift of galaxy A directly from stellar absorption lines will be a critical step to test the scenario.

Lastly, we note that the observed blueshift measured in the permitted broad lines should carefully be interpreted, because these have large uncertainties primarily due to broad absorption components and blending with other lines (Chiaberge et al. 2017). For example, such broad absorption is often associated with extreme outflow and seen in a number of local AGNs (Kaastra et al. 2014; Ebrero et al. 2016), making it challenging to attribute the observed blueshift solely to a recoiling event. An ideal line is Hβ, which is a good tracer of the kinematics of the accretion disk/inner broad-line region and is isolated from other emission lines; thus it can provide a clean measurement of the velocity shift of the accretion disk (Chiaberge et al. 2018). However, this line is located at ∼1 μm for the redshift (plus possible large blueshift) of 3C 186, and access to the entire line profile with current spectrographs may be challenging due to atmospheric absorption and the gap between CCD and IR detector. The IFU mode of NIRSpec on the James Webb Space Telescope (JWST) would be an ideal choice, allowing us to confirm the velocity shift and to locate spatial distribution simultaneously.

High-resolution submillimeter/millimeter observations of dense gas in 3C 186 would also be a critical step to confirm the physical origin of the blueshift. This is because dense gas in the molecular torus is considered to be orbiting at large distances and thus not gravitationally bounded to the recoiling black hole itself. Indeed, Decarli et al. (2014) tested this on a candidate, QSO J0927+2943 at z = 0.697, and identified molecular gas traced by the CO(2–1) line at the redshift of the quasar's broad line, arguing the need for another scenario to explain the properties of the quasar (see also Heckman et al. 2009; Shields et al. 2009).