A Trio of Massive Black Holes Caught in the Act of Merging

TL;DR

This paper reports the discovery of the first known triple active galactic nucleus system, combining multiwavelength observations to suggest these black holes may merge within 2 billion years, providing insights into galaxy evolution and gravitational wave sources.

Contribution

It presents the first confirmed triple Type 2 Seyfert nucleus with comprehensive multiwavelength data, advancing understanding of galaxy mergers and black hole evolution.

Findings

Identification of three active nuclei in a single galaxy system.

Multiwavelength data strongly support all three being active galactic nuclei.

Potential evolution into a bound triple black hole system within 2 Gyr.

Abstract

We report the discovery of SDSS J0849+1114 as the first known triple Type 2 Seyfert nucleus. It represents three active black holes that are identified from new spatially resolved optical slit spectroscopy using the Dual Imaging Spectrograph on the 3.5 m telescope at the Apache Point Observatory. We also present new complementary observations including the Hubble Space Telescope Wide Field Camera 3 U- and Y-band imaging, Chandra Advanced CCD Imaging Spectrometer S-array X-ray 0.5--8 keV imaging spectroscopy, and NSF Karl G. Jansky Very Large Array radio 9.0 GHz imaging in its most extended A configuration. These comprehensive multiwavelength observations, when combined together, strongly suggest that all three nuclei are active galactic nuclei. While they are now still at kiloparsec-scale separations, where the host-galaxy gravitational potential dominates, the black holes may evolve into a bound triple system in $\lesssim$2 Gyr. These triple merger systems may explain the overly massive stellar cores that have been observed in some elliptical galaxies such as M87, which are expected to be unique gravitational wave sources. Similar systems may be more common in the early universe, when galaxy mergers are thought to have been more frequent.

Tables

Table 1. Table 1: Host-galaxy Properties of SDSS J0849+1114

			$r$	$z$	$u-z$	$M_z$	$M_{z,c}$	$m_Y$	$m_U$	log$M_{\ast }$	log$M_{\ast ,\mathrm{ap}}$	log$M_{\bullet }$
ID	SDSS Designation	Redshift	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)
A	J084905.51+111447.2	0.0775	15.90	15.21	2.38	$-$22.52	$-$23.26	16.62	19.72	11.3	10.6	7.5
B	J084905.39+111445.8	0.0778	18.36	18.39	1.09	$-$19.35	$-$21.34	17.83	19.26	10.0	10.1	6.4
C	J084905.42+111451.0	0.0779	17.49	16.89	2.79	$-$20.26	$-$21.06	18.40	21.42	10.5	9.7	6.7

Table 2. Table 2: Emission-line properties of SDSS J0849+1114

							log $n_e$			${\sigma }_{\mathrm{gas}}$	${\sigma }_{\ast }$	log$M_{\bullet }$
ID	$\frac{\left[O\mathrm{III}\right]}{H\beta }$	$\frac{\left[N\mathrm{II}\right]}{H\alpha }$	$\frac{\left[OI\right]}{H\alpha }$	$\frac{\left[S\mathrm{II}\right]}{H\alpha }$	$\frac{H\alpha }{H\beta }$	$\frac{\left[S\mathrm{II}\right]\lambda 6717}{\left[S\mathrm{II}\right]\lambda 6731}$	(cm-3)	$\frac{\left[O\mathrm{III}\right]\lambda \lambda 4959,5007}{\left[N\mathrm{II}\right]\lambda \lambda 6548,6584}$	$\frac{\left[O\mathrm{III}\right]\lambda \lambda 4959,5007}{\left[S\mathrm{II}\right]\lambda \lambda 6717,6731}$	(km s-1)	(km s-1)	($M_{\odot }$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)
A†	7.92	0.592	0.200	0.445	4.75	1.15	2.49	6.87	13.06	134$\pm$3	193$\pm$19	7.9
B†	4.75	0.447	0.179	0.351	6.45	1.01	2.75	6.53	12.47	101$\pm$22	119$\pm$17	7.1
B‡	4.61	0.558	0.108	0.427	4.21	1.29	2.15	3.83	7.02	162$\pm$1	101$\pm$8	6.8
C†	4.78	0.493	0.140	0.442	5.00	1.38	1.77	4.66	7.50	113$\pm$14	144$\pm$53	7.4
C‡	2.42	0.600	0.108	0.479	4.89	1.05	2.68	2.14	3.84	162$\pm$1	109$\pm$17	6.9

Table 3. Table 3: HST Y 𝑌 Y -band Photometric Decomposition Results from GALFIT Analysis

	m${}_{Y,\mathrm{PSF}}$	m${}_{Y,S\acute{e}\mathrm{rsic}}$	Re
ID	(mag)	(mag)	(kpc)	n	B/T
(1)	(2)	(3)	(4)	(5)	(6)
A	18.4	15.7	2.4	1.3	0.076
B	18.9	16.1	8.2	1.5	0.071
C	21.4	17.0	3.9	3.5	-

Table 4. Table 4: X-ray measurements for SDSS J0849+1114

		Flux${}_{0.5-2\mathrm{keV}}$	Flux${}_{2-8\mathrm{keV}}$			$N_H$	log $L_{0.5-2\mathrm{keV}}$	log $L_{2-8\mathrm{keV}}$	log $L_{2-10\mathrm{keV}}$
ID	Counts	(10-6 photons cm-2 s-1)		HR	${\Gamma }_X$	(1022 cm-2)	(erg s-1)	(erg s-1)	(erg s-1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
A	178.8	6.76${}_{-0.73}{}^{}{}_{}{}^{+0.74}$	6.35${}_{-0.72}{}^{}{}_{}{}^{+0.73}$	$-$0.05${}_{-0.07}{}^{}{}_{}{}^{+0.07}$	1.12${}_{-0.48}{}^{}{}_{}{}^{+0.63}$	0.66${}_{-0.66}{}^{}{}_{}{}^{+1.18}$	41.67${}_{-0.04}{}^{}{}_{}{}^{+0.04}$	42.01${}_{-0.05}{}^{}{}_{}{}^{+0.04}$	42.13${}_{-0.05}{}^{}{}_{}{}^{+0.04}$
B	17.7	1.17${}_{-0.30}{}^{}{}_{}{}^{+0.35}$	0.77	$-$0.80${}_{-0.16}{}^{}{}_{}{}^{+0.06}$	-	-	41.27${}_{-0.13}{}^{}{}_{}{}^{+0.11}$	40.90	40.99
C	8.3	0.49${}_{-0.21}{}^{}{}_{}{}^{+0.26}$	0.76	$-$0.73${}_{-0.27}{}^{}{}_{}{}^{+0.06}$	-	-	40.90${}_{-0.23}{}^{}{}_{}{}^{+0.18}$	40.90	40.98

Table 5. Table 5: VLA radio measurements for SDSS J0849+1114

	S${}_{9.0\mathrm{GHz}}{}^{}{}_{}{}^{\mathrm{Peak}}$	S${}_{9.0\mathrm{GHz}}{}^{}{}_{}{}^{\mathrm{Int}}$	log L${}_{9.0\mathrm{GHz}}$	log L${}_{1.4\mathrm{GHz}}$	log L${}_{1.4\mathrm{GHz}}{}^{}{}_{}{}^{\mathrm{SF}}$	SFR${}_{1.4\mathrm{GHz}}$
ID	(mJy/beam)	(mJy)	(erg s-1 Hz-1)	(erg s-1 Hz-1)	(erg s-1 Hz-1)	($M_{\odot }$ yr-1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
A	2.58	5.52$\pm$0.17	30.09	30.49	$\lesssim$29.49	$\lesssim$20
B	0.018	0.018	27.61	28.08	28.08	0.8
C	0.348	0.472$\pm$0.014	29.03	29.43	$\lesssim$28.43	$\lesssim$2

Table 6. Table 6: Dust Extinction, Star Formation, and Intrinsic Luminosity Estimates

	$E\left(B-V\right)$	$A_u$	$A_z$	log $L_{u,c}$	SFR${}_{u,c}$		SFRd	log $L_{0.5-2\mathrm{keV}}^{\mathrm{SF}}$	log $L_{2-10\mathrm{keV}}^{\mathrm{SF}}$	log $L_{\mathrm{HX}}^{\mathrm{gal}}$	log $L_{\left[O\mathrm{III}\right],c}$
ID	(mag)	(mag)	(mag)	($\frac{\mathrm{erg}}{s\mathrm{Hz}}$)	($\frac{M_{\odot }}{\mathrm{yr}}$)	$D$(4000)	($\frac{M_{\odot }}{\mathrm{yr}}$)	(erg s-1)	(erg s-1)	(erg s-1)	(erg s-1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
A†	0.51	2.50	0.74	29.0	8	-	-	40.57	40.61	40.49	43.10
B†	0.82	4.01	1.19	29.8	70	-	-	41.48	41.53	41.06	42.72
B‡	0.39	1.90	0.56	29.0	7	1.14	0.2	[39.0, 40.5]	[39.0, 40.5]	[39.1, 40.1]	41.66
C†	0.56	2.75	0.81	28.5	2	-	-	39.88	39.95	39.79	42.39
C‡	0.54	2.64	0.78	28.4	2	1.24	0.1	[38.7, 40.0]	[38.7, 40.0]	[39.5, 39.8]	41.63

Equations

$$\Sigma(r)=\Sigma_{e}\,{\rm exp}\bigg{[}-\kappa\bigg{(}\Big{(}\frac{r}{r_{e}}\Big{)}^{1/n}-1\bigg{)}\bigg{]}$$

$$\mathrm{lo}{\mathrm{g}}_{10}\left(\frac{M_{\ast }}{M_{\odot }}\right)=-0.4(M_z-M_{z,\odot })-0.179+0.151(u-z),$$

$$\frac{{\rm SFR}_{u,\,{\rm c}}}{M_{\odot}~{}{\rm yr}^{-1}}=\bigg{(}\frac{L_{u,\,{\rm c}}}{1.81\times 10^{28}~{}{\rm ergs~{}s}^{-1}{\rm Hz}^{-1}}\bigg{)}^{1.186},$$

$$\begin{array}{r}{L}_{0.5-2\mathrm{keV}}^{\mathrm{SF}}=4.5\times 10^{39}\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}{\mathrm{ergs}\mathrm{s}}^{-1},\\ L_{2-10\mathrm{keV}}^{\mathrm{SF}}=5.0\times 10^{39}\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}{\mathrm{ergs}\mathrm{s}}^{-1},\end{array}$$

$$L_{\mathrm{HX}}^{\mathrm{gal}}=\alpha M_{\ast }+\beta \mathrm{SFR},$$

$$U_i\equiv \frac{Q}{4\pi d_i^2{n}_{e}c},$$

$$\mathrm{log_{10}}\bigg{(}\frac{M_{\bullet}}{M_{\odot}}\bigg{)}=(7.45\pm 0.08)+(1.05\pm 0.11)\,\mathrm{log_{10}}\bigg{(}\frac{M_{\ast}}{10^{11}M_{\odot}}\bigg{)},$$

$$\mathrm{log_{10}}\bigg{(}\frac{M_{\bullet}}{M_{\odot}}\bigg{)}=(7.96\pm 0.03)+4.02\,\mathrm{log_{10}}\bigg{(}\frac{\sigma_{\ast}}{200\,{\rm km}\,{\rm s}^{-1}}\bigg{)},$$

$$t^{{\rm gal}}_{{\rm fric}}=\frac{2.7}{\ln\Lambda}\frac{r}{30\,{\rm kpc}}\bigg{(}\frac{\sigma_{{\rm H}}}{200\,{\rm km\,s^{-1}}}\bigg{)}^{2}\bigg{(}\frac{100\,{\rm km\,s^{-1}}}{\sigma_{{\rm S}}}\bigg{)}^{3}~{}{\rm Gyr},$$

$$t^{{\rm BH}}_{{\rm fric}}=\frac{19}{\ln\Lambda}\bigg{(}\frac{r}{5\,{\rm kpc}}\bigg{)}^{2}\frac{\sigma_{{\rm H}}}{200\,{\rm km\,s^{-1}}}\frac{10^{8}\,M_{\odot}}{M_{\bullet}}~{}{\rm Gyr},$$

Full text

Galaxy mergers – Galaxy triplets – Active galaxies – AGN host galaxies – X-ray active galactic nuclei — High energy astrophysics – Radio cores – Black hole physics – Astrophysical black holes – Surveys

††software: CASA (McMullin et al., 2007), CIAO (Fruscione et al., 2006), Sherpa (Freeman et al., 2001; Doe et al., 2007), ChiPS (Germain et al., 2006), pPXF (Cappellari & Emsellem, 2004), qsofit (Guo et al., 2018; Shen et al., 2019)

1 Introduction

Most local galaxies, if not all, are believed to harbor a massive black hole (BH) at the center222There are, however, galactic nuclei that do not appear to contain massive black holes (MBHs), such as M33 (Gebhardt et al., 2001). (Kormendy & Ho, 2013). Black hole pairs and triples are natural outcomes of the hierarchical galaxy and BH assembly process. Studies based on the Sloan Digital Sky Survey (SDSS; York et al., 2000) galaxies show that ${\lesssim}2$% of the (major) merging systems involve three galaxies (e.g., Darg et al., 2010). While hundreds of active BH pairs (i.e., dual active galactic nuclei) and candidates have been found in galaxy mergers (e.g., Liu et al., 2011b), only three kiloparsec-scale (kpc-scale for short) triple-AGN candidates are known (Barth et al., 2008; Liu et al., 2011a; Schawinski et al., 2011). However, their identification based on optical diagnostic emission-line ratios has been inconclusive. Alternative ionizing sources remain possible, such as stellar and shock-heating processes.

The galaxy SDSS J084905.51$+$111447.2 (hereafter SDSS J0849+1114 for short) at redshift $z=0.078$ contains three optical stellar nuclei within a projected $\sim$5 kpc radius that show a disturbed morphology (Figure 1). SDSS J0849+1114 was identified as a candidate to host a kpc-scale triple MBH from the largest sample of optically selected candidate AGN pairs (Liu et al., 2011b) based on the SDSS Seventh Data Release (Abazajian et al., 2009). The SDSS has two spectra taken with a 3*′′* diameter fiber centered on Nucleus B (projected at 2$\farcs$3 or 3.4 kpc to the southwest of A) and C (projected at 3$\farcs$6 or 5.2 kpc to the northwest of A), classifying B as a Type 2 Seyfert and C as a low-ionization nuclear emission-line region333C is classified as a LINER by Liu et al. (2011b) based on emission-line ratios measured from the SDSS spectrum because the Ho et al. (1997) criterion is adopted to separate Seyferts from LINERs; it would have been classified as a Seyfert if the Schawinski et al. (2007) criterion were adopted. Our new APO/DIS measurements unambiguously classify C as a Seyfert regardless of which criterion is used. (LINER; Liu et al., 2011b).

To understand the physical nature of Nucleus A and to more clearly separate light coming from each individual nucleus, we here present spatially resolved follow-up optical spectroscopy with the Dual Imaging Spectrograph (DIS) on the 3.5m telescope at the Apache Point Observatory. After carefully separating the light from each nucleus and subtracting the host-galaxy stellar continuum (Figure 2), all three nuclei are optically classified as Type 2 (obscured) Seyferts based on the diagnostic narrow emission-line ratios according to the classical BPT diagram (Figure 3).

SDSS J0849+1114 represents the first known case of a triple Type 2 Seyfert nucleus. Unlike SDSS J0849+1114, previous triple-AGN candidates all contain at least one nucleus classified as either LINER or AGN-H II composite, which may be due to excitation by evolved stars, star formation, and/or shock heating instead of AGN ionization. For example, the disk galaxy NGC 3341 contains three nuclei that are optically classified as a Seyfert, a LINER, and a LINER or LINER-H II composite (Barth et al., 2008). X-ray and radio follow-up observations suggest that both the primary galaxy and one of the secondary dwarf companions are unlikely to host an AGN (Bianchi et al., 2013). Another potential candidate was reported at $z=1.35$ based on Hubble Space Telescope (HST) grism spectroscopy (Schawinski et al., 2011), although the spectral resolution was too low to reliably resolve the individual components, which may be from a clumpy star-forming disk galaxy rather than a triple AGN. Yet another candidate, now disputed, J1502+1115 (at $z=0.39$), was identified by the Very Long Baseline Interferometry (VLBI; Deane et al., 2014), although its two nuclei were later found to be radio hot spots instead of a pair of BHs (Wrobel et al., 2014b). Finally, in another candidate from the SDSS candidate AGN pair sample, SDSSJ1027+1749, only one of the nuclei is optically classified as a Seyfert, whereas the other two are a LINER and a composite (Liu et al., 2011a).

While the triple Type 2 Seyfert nucleus suggests three active MBHs in SDSS J0849+1114, optical diagnostics may be inconclusive. There could be only one or two active BHs ionizing all three merging components, resulting in three Seyfert-like nuclei, although this scenario is unlikely considering their positions on the BPT diagram. To further determine the excitation mechanism of the triple Type 2 Seyfert nucleus in SDSS J0849+1114, here we also present a comprehensive follow-up campaign including HST Wide Field Camera 3 (WFC3) $U$- and $Y$-band imaging, Chandra Advanced CCD Imaging Spectrometer S-array (ACIS-S) X-ray 0.5–8 keV imaging spectroscopy, and NSF Karl G. Jansky Very Large Array (VLA) radio 9.0 GHz imaging in its most extended A configuration.

The paper is organized as follows. §2 presents our target selection and new observations from APO 3.5m/DIS (§§2.2, 2.3), HST/WFC3 (§§2.4, 2.5), Chandra ACIS-S (§2.6), and VLA A-config (§2.7). §3 shows our results on the host-galaxy internal dust extinction (§3.1), stellar mass (§3.2) and nuclear star formation rate (SFR; §3.3), X-ray contribution from nuclear star formation (§3.4), photoionization diagnostics based on optical spectra (§3.5), and the nature of the nuclear ionizing sources (§3.6). We then discuss the implications of our results in §4. Finally, we summarize our main findings and discuss future work in §5. Throughout this paper, we assume a concordance cosmology with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and $H_{0}=70$ km s*-1* Mpc*-1* and use the AB magnitude system (Oke, 1974).

2 Observations, Data Reduction, and Data Analysis

2.1 Target Selection

SDSS J0849+1114 (Table 1) was selected from an [O III]-selected AGN pair sample (Liu et al., 2011b) systematically identified in the SDSS DR7 (Abazajian et al., 2009). We searched for kpc-scale triple-AGN candidates among kpc-scale AGN pairs with a companion galaxy within 10 kpc in projection from the first two nuclei. We then conducted spatially resolved optical spectroscopy to measure redshifts and to estimate excitation mechanisms for all nuclei. While in most systems the third nucleus turned out to be an inactive neighbor, we found three candidate kpc-scale triple AGNs (Liu et al., 2011a) for further study. SDSS J0849+1114 is unique among the three candidates in that it is the only one with a triple Type 2 Seyfert nucleus. SDSS J0849+1114 was also independently identified as a dual AGN/triple-AGN candidate (Pfeifle et al., 2019a) in a sample of IR-selected mergers based on the Wide-field Infrared Survey Explorer data (Wright et al., 2010). After we submitted this manuscript, another work (Pfeifle et al., 2019b) also independently proposed SDSS J0849+1114 as a triple AGN, although we argue that the comprehensive evidence shown in the present paper (in particular the stringent constraints on the alternative scenario of a nuclear starburst based on the high-resolution VLA and HST $U$-band imaging) is significantly stronger.

We distinguish between kpc-scale triple AGN (Liu et al., 2011a, i.e., the focus of this work) from physical triple quasars (Djorgovski et al., 2007; Farina et al., 2013; Hennawi et al., 2015). The former have typical projected physical separations of a few kiloparsec and velocity offsets of $<$300 km s*-1* when the host galaxies are already in direct gravitational interactions. The latter have typical projected physical separations of a few hundred kiloparsec and velocity offsets of $<$1000 km s*-1*, suggesting that the quasars belong to the same dark matter halo, although it is unclear if their host galaxies are already in direct gravitational interaction with each other.

2.2 APO 3.5m/DIS Long-slit Spectroscopy

We obtained long-slit spectra for SDSS J0849+1114 on the night of 2011 March 6 UT using the Dual Imaging Spectrograph (DIS) on the Apache Point Observatory 3.5 m telescope. The sky was non-photometric, and the seeing ranged from 0$\farcs$8 to 2$\farcs$0, with a median of $\sim$1$\farcs$2. The field of view of DIS was 4${}^{\prime}\times$6*′* with a pixel size of 0$\farcs$414. We used a 1$\farcs$5$\times$6*′* slit with the B1200+R1200 gratings centered at 510 and 700 nm. The spectral coverage was 450–560 nm and 640–760nm with an instrumental resolution of $\sigma_{{\rm inst}}\approx$60 and 30 km s*-1* and a dispersion of 0.62 and 0.58 Å pixel*-1* in the blue and red channels. The slit was oriented at a position angle of 43*∘* and 159*∘* to go through Nucleus A and B, and A and C (Figure 6). The total effective exposure time was 5400 s, 2700 s, and 2700 s for A, B, and C. We observed standard stars G191B2B and HZ44 to calibrate spectrophotometry.

We reduced the DIS data following standard IRAF444IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation. procedures (Tody, 1986). We extracted 1D spectra using a 2*′′* diameter aperture for each nucleus. We applied a telluric correction from standard stars to the extracted 1D spectra. The median signal-to-noise ratio (S/N) we achieved is $\sim$10–40 pixel*-1* (Figure 2).

2.3 Optical Spectral Analysis

We first measure the host-galaxy stellar continuum in the emission-line free regions of the observed spectrum and then model the emission-line flux over the host-subtracted spectrum. For the host-galaxy spectral fitting, we use the penalized pixel-fitting (pPXF) method (Cappellari & Emsellem, 2004), which works directly in the pixel space using the maximum penalized likelihood formalism to extract the most information from the spectra while suppressing noise in the solution. For each nucleus, we measured the host redshift by fitting the continuum with galaxy templates produced by population synthesis models (Bruzual & Charlot, 2003) following the procedure of Liu et al. (2009). We convolved the galaxy templates with the stellar velocity dispersion $\sigma_{\ast}$ measured over small spectral windows containing strong stellar absorption features.

After subtracting the stellar continuum using the best-fit pPXF model, we fit the narrow emission lines simultaneously with multiple Gaussian models constrained to have the same velocity and line width using the spectral fitting code qsofit (Guo et al., 2018; Shen et al., 2019). We also verified the host redshift by measuring the velocity of the combined set of narrow emission lines.

Figure 2 shows the rest-frame 1D spectrum and our best-fit models for the stellar continuum and the continuum-subtracted emission lines. Table 2 lists the emission-line measurements from our best-fit models.

2.4 HST WFC3 $U$- and $Y$-band Imaging

We observed SDSS J0849+1114 using the WFC3 on board the HST on 2013 March 18 UT (Program GO 13112; PI Liu). Images were taken in the UVIS/F336W ($U$ band, with pivot $\lambda_{{\rm p}}=$335.5 nm and width of 51.1 nm; Dressel 2010) and IR/F105W (wide $Y$ band, with $\lambda_{{\rm p}}=$1055.2 nm and width of 265.0 nm, Dressel 2010) filters in a single orbit. The total net exposure times were 2151 s in $U$ and 239 s in $Y$. To cover SDSS J0849+1114 and the nearby field adequately for background subtraction, we adopted a 1k$\times$1k subarray in $U$ for a field of view (FOV) of 40${}^{\prime\prime}\times$37*′′* and a 512$\times$512 subarray in $Y$ for an FOV of 72${}^{\prime\prime}\times$64*′′*. The WFC3 UVIS CCD has a pixel size of 0$\farcs$039, and the IR detector has a pixel size of 0$\farcs$13. We dithered the observations to properly sample the point-spread functions (PSFs) while rejecting cosmic rays and bad pixels.

We reduced the HST data following standard procedures. Geometric distortion and pixel area effects were corrected for. We combined dithered frames while rejecting cosmic rays and hot pixels. The final image has a pixel scale of 0$\farcs$06 for the $Y$-band image to be Nyquist-sampled. Typical relative astrometric accuracy is 0$\farcs$004 for UVIS and 0$\farcs$01 for IR images. The absolute astrometric accuracy of the reduced WFC3 image is limited by the positional uncertainty of guide stars and the calibration uncertainty of the fine guidance sensor to the instrument aperture. For better absolute astrometric accuracy, we registered the $Y$-band image to the SDSS astrometry (Pier et al., 2003) using reference objects selected from the FOV with measured SDSS astrometry. The absolute astrometric uncertainty of the registered image is estimated as $\sim$0$\farcs$15, which includes $\sim$0$\farcs$055, the typical SDSS astrometric uncertainty (Pier et al., 2003), and $\sim$0$\farcs$085, the statistical uncertainty of the standard deviation of the WCS fit. We tie the $U$-band astrometry to the $Y$ band because of the small relative astrometric uncertainty between adjacent exposures. Figure 1 shows the reduced and calibrated HST images for SDSS J0849+1114.

2.5 Analysis of Host-galaxy Morphology and Structural Decomposition

We use the HST $Y$-band image to study host-galaxy morphology and old stellar populations. Given the highly disturbed morphology in the complex merging galaxy, we use GALFIT (Peng et al., 2010) to model the 2D surface brightness profile with multiple structural components. Instead of achieving a perfect fit to the data, our goal is to try and decompose the three galaxies from each other and to try and decompose the disk and bulge components, if any, in each galaxy. We use stars within the FOV to model the PSF and employ a constant sky background in the model. For each galaxy, we use one S${\rm\acute{e}}$rsic model for a disk or a resolved bulge component and a PSF model for any unresolved component (e.g., from a compact stellar bulge, a nuclear starburst, or continuum emission from an unresolved narrow line region). The S${\rm\acute{e}}$rsic model profile is given by

[TABLE]

where $\Sigma(r)$ is the pixel surface brightness at radial distance $r$, $\Sigma_{e}$ is the pixel surface brightness at the effective radius $r_{e}$, and $\kappa$ is a parameter related to the S${\rm\acute{e}}$rsic index $n$. $n=1$ for an exponential profile, while $n=4$ (corresponding to $\kappa=7.67$) for a de Vaucouleurs profile. Bulge-dominated galaxies have high $n$ values (e.g., $n>2$), while disk-dominated galaxies have $n$ close to unity. We also tried using an additional S${\rm\acute{e}}$rsic for any resolved bulge, but adding another component did not significantly improve the fit.

Figure 4 shows our best-fit GALFIT model. The residual image still shows significant patterns, which is expected in highly disturbed mergers. Figure 5 shows the radial surface brightness profiles for individual components. Table 3 lists the best-fit model parameters. Our best-fit model for Nucleus A is dominated by a disk component with S${\rm\acute{e}}$rsic index $n=1.3$. Nucleus B is dominated by an extended disk component with S${\rm\acute{e}}$rsic index $n=1.5$, which is likely being tidally striped by A. Nucleus C is dominated by a bulge component with S${\rm\acute{e}}$rsic index $n=3.5$, although it does show an extended disk component in the residual image, which is again likely being tidally striped by A.

2.6 Chandra ACIS-S X-Ray 0.5–8 keV Imaging Spectroscopy, X-Ray Data Reduction, and Data Analysis

We observed SDSS J0849+1114 with the ACIS-S on board the Chandra X-ray Observatory on 2013 February 8 UT (Program GO3-14104A, ObsID=14969; PI Liu). SDSS J0849+1114 was observed on-axis on the S3 chip with an exposure time of 20 ks. Following our Cycle 14 observation, there was a Cycle 17 ACIS-S observation, taken on 2016 March 3 UT with an exposure time of 21 ks (ObsID=18196; PI Satyapal). We combine both observations to maximize the S/N. The total effective exposure time of the combined image is 41 ks.

We reprocessed the Chandra data using CIAO v4.8 with the calibration files CALDB v4.7.0 following standard procedures. The light curve of each observation was examined to ensure that there was no significant particle background. We generate counts maps, exposure maps and PSF maps at the original pixel scale in all bands (0$\farcs$492 pixel*-1*). We weigh exposure maps with a fiducial model for Type 2 Seyferts (Green et al., 2009) assuming an absorbed power law with a photon-index $\Gamma_{{\rm X}}$=1.7 and a column density $N_{{\rm H}}=$1022 cm*-2*. The CIAO tool wavdetect is used to search for tentative X-ray sources coincident with the three nuclei. For each resultant source, we derive a background-subtracted and exposure-map–corrected photon flux from a 1*′′* radius circle (representing $\sim$ 90% of the enclosed energy) using the CIAO tool aprate. This extraction aperture ensures that there is no overlap between different nuclei. The local background is extracted typically between 2*′′* and 5*′′*, and other detected point sources in this region are removed. Nuclei with 3$\sigma$ photon flux lower limits greater than 0 are considered as a firm detection. For the non-detected nuclei, we estimate their 3$\sigma$ photon flux upper limits. Nucleus A is detected in both soft and hard bands, whereas B and C are detected in the soft band only. Figure 6 shows the ACIS image of SDSS J0849+1114 in the soft (0.5–2 keV) and hard (2–8 keV) bands. Given the low count level, we do not apply any smoothing to avoid artifacts. All X-ray detections are consistent with the HST $Y$-band nucleus positions within astrometric uncertainties.

For Nucleus A, which has sufficient counts (76 and 82 net counts in two observations), we perform a spectral analysis, adopting C-statistics, which is suitable for the low-count regime (Cash, 1979). No significant spectral variation is found between the two observations, hence the two spectra are jointly fitted with all model parameters linked. We first try a single absorbed power-law model, i.e., $\rm{phabs\times(zphabs\times pow)}$ in XSPEC, where $\rm{phabs}$ represents the Galactic foreground absorption (fixed at 3$\times$1020 cm*-2*) and $\rm{zphabs}$ accounts for internal absorption to the AGN. This model, however, provides a poor fit to the spectra. Thus we introduce an additional power-law component to account for potential contribution from unresolved X-ray binaries and/or hot gas within the extraction aperture, i.e., $\rm{phabs\times(zphabs\times pow+pow)}$ in XSPEC. This phenomenological model provides an acceptable fit to the spectra, and the best-fit photon-index and absorption column density are typical of Type 2 Seyferts. Figure 7 shows the spectra and best-fit model for A. Table 4 lists the X-ray measurements.

For Nuclei B and C, which have too few net counts for a spectral analysis, we estimate their X-ray spectral properties using hardness ratios (HR), defined as ${\rm HR}\equiv(H-S)/(H+S)$, where $H$ and $S$ are the number of net counts in the hard (2–8 keV) and soft (0.5–2 keV) X-ray bands. We use the Bayesian estimation of hardness ratios (Park et al., 2006) to estimate HRs and uncertainties, appropriate for the low-count regime. Table 4 lists the HR results and the estimated X-ray luminosity for B and C. We assume an absorbed power law with a photon index of $\Gamma_{{\rm X}}=1.7$ and an intrinsic absorption column density of $N_{\rm H}$=1022 cm*-2* to calculate the luminosities for B and C based on their net count rates. The measured HR for B and C suggest softer spectra than the model assumption, although the inferred luminosity difference is insignificant. Strictly speaking, the estimated X-ray luminosities should be taken as lower limits considering the possibility of significantly underestimating the true absorbing column in the dusty merger due to the poor constraints from the low X-ray counts.

2.7 VLA Radio 9.0 GHz Imaging, Radio Data Reduction, and Data Analysis

SDSS J0849+1114 has an integrated 1.4 GHz flux density of 38.6$\pm$1.2 mJy from NVSS (Condon et al., 1998) and 35.4$\pm$1.8 mJy555FIRST error bars need to be increased by 5% in quadrature (Perley & Butler, 2017). from FIRST (White & Becker, 1992). The FIRST source is centered on Nucleus A with a peak flux density of 30.2$\pm$1.5 mJy/beam, but the resolution is insufficient to separate emission from each individual nucleus.

To resolve the triple nucleus, we observed SDSS J0849+1114 with the VLA in its A configuration at 9.0 GHz on 2012 October 6 under project code SD0279 (PI Liu). The total exposure time was 1 hr with an on-source time of 0.5 hr. A total of 27 antennas were used. We reduced the new VLA data following standard procedures (Wrobel et al., 2014a) using the Common Astronomy Software Aplications (CASA) package (McMullin et al., 2007). The beam size is 0$\farcs$19$\times$0$\farcs$19. Nuclei A and C were detected with peak flux densities of 2.58 mJy/beam and 0.348 mJy/beam and integrated flux densities of 5.52$\pm$0.17 mJy and 0.472$\pm$0.014 mJy, respectively.666VLA error bars need to be increased by 3% in quadrature (Perley & Butler, 2017). We have also detected some features with peak flux densities of 0.05–0.07 mJy/beam between Nucleus A and C, which could be galactic-scale star formation and/or AGN outflows. A 3$\sigma$ upper limit of 0.018 mJy/beam was set for Nucleus B.

We extrapolate our VLA 9.0 GHz measurements to 1.4 GHz assuming $f_{\nu}\propto\nu^{-0.5}$, which is canonical for Seyfert galaxies (Ho & Peng, 2001). Our approach is conservative because synchrotron self-absorption could produce a flatter spectral index. A flatter spectral index would result in a lower 1.4 GHz flux. This would in turn translate into a lower radio-based SFR estimate, which would make the case for an AGN even stronger (i.e., to explain the observed X-ray emission). For Nucleus B, we also explore an alternative model assuming $f_{\nu}\propto\nu^{-0.59}$ which is appropriate for star-forming galaxies (Klein et al., 2018), assuming the emission is dominated by nonthermal synchrotron. Again our approach is conservative because we assume the average spectral index for the nonthermal synchrotron emission observed in star-forming galaxies. Adopting the thermal free-free emission instead would result in a much flatter spectral index. Table 5 lists our radio measurements.

3 Results

3.1 Host-galaxy Internal Dust Extinction

We estimate the extinction correction using the Balmer decrement method (Osterbrock, 1989) based on the emission-line intensity ratio H$\alpha$/H$\beta$ measured from our spatially resolved APO/DIS spectra (Table 2). We assume a dust-screen model with the intrinsic Case B value of H$\alpha$/H$\beta$=2.87 for $T=10^{4}$ K and the extinction curve of Cardelli et al. (1989) with $R_{V}$=3.1. To robustly measure H$\alpha$/H$\beta$, we have carefully subtracted the host-galaxy stellar continuum to minimize its error induced by Balmer absorption (Figure 2).

Table 6 lists the estimated color excess and SDSS $u$- and $z$-band dust extinctions for each nucleus. Typical 1$\sigma$ statistical errors ($<$0.1 mag) are insignificant compared to systematic errors, which are difficult to quantify given the unknown dust geometry. Our baseline dust-screen model is simplistic. We neglect the effects of metallicity and ionization state in the approximation of assuming a fixed unattenuated Case B ratio. The metallicity-dependence of the Case B ratio would cause an overestimate of the dust attenuation by up to $\sim$0.5 mag at rest frame 6563 Å for the most metal rich galaxies (Brinchmann et al., 2004).

Table 6 also lists the measurements based on the SDSS fiber spectra for Nuclei B and C. While the SDSS spectroscopic aperture is larger than that our APO/DIS spectra (and is therefore more contaminated by circumnucleus emission), the SDSS wavelength coverage is larger than that of APO/DIS spectra. In particular, the SDSS spectrum of Nucleus B shows strong Balmer absorption from a post-starburst stellar population (Figure 8). The strong Balmer absorption series in the blue end of the spectrum are crucial in helping to break the degeneracy between Balmer absorption and Balmer emission. On the other hand, our DIS spectrum does not cover the blue side, which contains the Balmer absorption series. Because of the limited coverage, it is likely that we have underestimated the absorption components (and overestimated the emission components) in H$\alpha$, H$\beta$, and H$\gamma$ in the DIS measurements. While the effects on the BPT ratios are likely to be minor, this would bias the Balmer decrement (Table 2), and by extension, the dust extinction (Table 6) estimate high, resulting in an overestimated SFR estimate. In view of this caveat, we adopt the SDSS-based measurements as our fiducial estimates for Nucleus B.

3.2 Host-galaxy Stellar Mass

We estimate the total stellar mass for the host of each nucleus using the $k$-corrected SDSS $z$-band absolute magnitude at rest frame after correction for dust extinction, $M_{z,\,{\rm c}}$, and the stellar mass-to-light ratio $M/L$ inferred from SDSS $u-z$ color ($k$-corrected to the rest frame). We adopt the empirical relation given by Bell et al. (2003)

[TABLE]

where $M_{*}$ is the galaxy stellar mass in solar units and $M_{z,\odot}=4.51$ is the absolute magnitude of the Sun in the $z$ band (Blanton & Roweis, 2007). Typical systematic uncertainties from galaxy age, dust, and bursts of star formation are $\sim$0.1 dex for optical $M/L$ and are larger for bluer galaxies (Bell et al., 2003). Uncertainties in stellar evolution are likely to dominate the systematic error in stellar mass, which is expected to be $\sim$0.3 dex (Conroy, 2013).

We also estimate the stellar mass contained within the same region as the X-ray extraction aperture for each nucleus, $M_{\ast,\,{\rm ap}}$. Along with the specific SFR inferred from optical continuum spectral index, $M_{\ast,\,{\rm ap}}$ is used to obtain an independent estimate for the SFR (see below for details). We estimate $M_{\ast,\,{\rm ap}}$ using a similar approach for estimating the total stellar mass, except that we use the aperture magnitudes $m_{U}$ and $m_{Y}$ and the aperture $u-z$ color estimated by applying $k$-corrections to $m_{U}$ and $m_{Y}$. Assuming a flat local spectrum between $u$ and $U$ and between $z$ and $Y$, the $k$-corrections are given by $m_{z}$ = $m_{Y}$ + 0.36 and $m_{u}$=$m_{U}$ $-$ 0.13. Table 1 lists our results.

3.3 Host-galaxy Nuclear SFR

We estimate the host-galaxy nuclear SFR from the HST-measured $U$-band luminosity following the method of Liu et al. (2013). We choose the HST $U$-band over the SDSS $u$-band imaging because the former resolves the emission from each nucleus better. We measure the HST $U$-band magnitude for each nucleus using aperture photometry. We adopt the same aperture as in the X-ray source extraction to minimize any aperture effects. We then convert the HST $U$-band magnitude to rest-frame SDSS $u$-band using a $k$-correction estimated by comparing the SDSS $u$-band aperture magnitude with the HST $U$-band magnitude from the same aperture as the SDSS. Table 6 lists $L_{u,\,{\rm c}}$, the intrinsic $u$-band luminosity density corrected for internal dust extinction, and the inferred SFR.

The $U$ band covers the rest frame 2870–3380 Å for SDSS J0849+1114. For Type 2 Seyferts where the AGN power law is obscured in the UV/optical, the luminosity in this range is dominated by continuum photosphere emission of young stellar populations in the host galaxy, providing a useful indicator for SFR (Cram et al., 1998). It is relatively free from strong emission lines from ionized gas, AGN light from the obscured nucleus scattered into our line of sight by dust and/or gas (Zakamska et al., 2006), considering the moderate AGN luminosity and the absence of a broad H$\beta$ or a broad H$\alpha$ in the optical spectrum (Liu et al., 2009), and the nebular continuum from ionized gas in the AGN emission-line region (Osterbrock, 1989). Accounting for possible contamination from all of these non-star formation related processes would lower the $U$-band–based SFR estimates and strengthen the case of AGN excitation.

To infer the SFR from $L_{u,\,{\rm c}}$, we have adopted the empirical calibration (Hopkins et al., 2003) obtained based on 2625 SDSS star-forming galaxies, given by

[TABLE]

which is valid for 2$\times$ 1028 ergs s*-1* Hz*-1* ${\lesssim}L_{u,\,{\rm c}}{\lesssim}$1030 ergs s*-1* Hz*-1*. Equation 3 has been calibrated against H$\alpha$-based SFR with an rms scatter of 0.13 dex. The calibrated $u$-band SFR has been shown to be consistent with the radio 1.4 GHz based SFR with an rms scatter of 0.23 dex (Hopkins et al., 2003).

We compare the $U$-band-based SFR with an independent estimate inferred from the 4000 Å break $D$(4000) using the narrow definition (Balogh et al., 1998). SFR indicators based on optical emission lines such as H$\alpha$ luminosity cannot be easily applied to SDSS J0849+1114 given the likely significant contamination from AGN excitation. We adopt $D$(4000) because all three nuclei in SDSS J0849+1114 are Type 2, i.e., obscured Seyferts, in which the AGN power law is obscured in the optical and the observed continuum is dominated by host-galaxy starlight. From $D$(4000), we infer the specific SFR, i.e., SFR/$M_{\ast}$, assuming the empirical calibration based on SDSS star-forming galaxies (Brinchmann et al., 2004). We then estimate the SFR by multiplying the specific SFR with the stellar mass within the same region as the X-ray extraction aperture, $M_{\ast,\,{\rm ap}}$ (Table 1). The $D$(4000) based SFR, SFRd, is consistent with the $U$-band-based SFR within uncertainties (Table 6), lending support to our SFR estimates.

For Nucleus B, we adopt the SDSS-based SFR estimates as our fiducial value. The DIS-based SFR is likely to be significantly overestimated due to the likely overestimated dust extinction correction. The more moderate SFR estimate is also more in line with the total (i.e., from three nuclei combined) SFR (${\sim}$20 M⊙ yr*-1*) inferred from the total IR luminosity (log($L_{{\rm IR}}/L_{\odot}$=11.43$\pm$0.03; Pfeifle et al., 2019a) as well as the strict upper limit from its radio luminosity density (see §2.7 for details), both of which are much less subject to uncertainties caused by the dust-attenuation correction.

3.4 X-Ray Contribution from Nuclear Star Formation

The estimated unabsorbed hard X-ray luminosities for the three nuclei are close to or below 1042 erg s*-1*, the characteristic upper limit for X-ray luminous starburst galaxies (Zezas et al., 2001). To determine whether the X-ray luminosity is dominated by star-formation-induced processes or by AGN excitation, X-ray spectral properties provide a discriminating diagnostic, although uncertainties in the X-ray measurements are too large to draw a firm conclusion.

Alternatively, we derive independent estimates on the expected X-ray luminosity due to star-formation-related processes in each nucleus and compare with the total X-ray luminosity to assess whether AGN is needed as an additional excitation source. We estimate the X-ray contribution from nuclear star formation using SFR${}_{u,\,{\rm c}}$, the $U$-band-based aperture SFR corrected for dust extinction. We adopt the empirical relation given by

[TABLE]

which is calibrated using 23 nearby star-forming galaxies (Grimm et al., 2003; Ranalli et al., 2003). The rms scatter is 0.27 dex in the soft and 0.29 dex in the hard X-ray band. Table 6 lists the resulting estimates.

Our estimates of the expected X-ray luminosity due to star-formation-related processes in the nucleus may be highly uncertain due to systematics in the model assumptions. The uncertainty is likely dominated by systematic errors in the $U$-band extinction correction. The Balmer decrement is estimated using emission-line measurements for each nucleus from our APO/DIS spectra. Given the lower angular resolution and projection effects, the aperture of the emission-line measurement does not exactly match with the HST $U$-band measurement, causing errors due to aperture mismatch. The reddening uncertainty due to the aperture mismatch is likely to be $<$1 mag (Liu et al., 2013). More importantly, however, the color excesses we derived using the Balmer decrement method may not represent the true dust attenuation. For example, dust could concentrate on scales smaller than where the Balmer lines are emitted. Patchy and optically thick dust clouds could result in high extinction with little reddening.

Our inferred star-formation-related X-ray luminosity is not sensitive to the initial mass function (IMF) uncertainty, assuming no large systematic IMF variations among star-forming galaxies (Scalo, 1986; Kroupa, 2001; Chabrier, 2003). A Salpeter IMF (Salpeter, 1955) is assumed with mass limits of 0.1 $M_{\odot}$ and 100 $M_{\odot}$ in the calibration of SFR from $U$-band luminosity (Equation 3). The adopted SFR–X-ray-luminosity relation (Equation 4; Ranalli et al., 2003) is calibrated under the same assumptions about the IMF and mass range (Kennicutt, 1998).

We double-checked our results using the empirical correlation between galaxy-wide 2–10 keV luminosity L${}^{{\rm gal}}_{{\rm HX}}$ (assumed to be the combined emission from high-mass X-ray binaries and low-mass X-ray binaries), SFR, and stellar mass from Lehmer et al. (2010). This calibration is based on a sample of 17 local luminous IR galaxies. The relation is given by

[TABLE]

where $\alpha=(9.05\pm 0.37)\times 10^{28}$ erg s*-1* $M_{\odot}^{-1}$ and $\beta=(1.62\pm 0.22)\times 10^{39}$ erg s*-1* ($M_{\odot}$ yr*-1*)-1. The stellar mass and SFR-derived hard X-ray luminosities are consistent with the SFR-derived hard X-ray estimates within the uncertainties (Table 6), lending confidence to our estimates.

3.5 Photoionization Estimates

We use [O III] $\lambda\lambda$4959,5007/[N II] $\lambda\lambda$6548,6584 and [O III] $\lambda\lambda$4959,5007/[S II] $\lambda\lambda$6717,6731 (corrected for reddening) to quantify the level of ionization in the emission-line gas. Both line flux ratios increase as the level of ionization increases. The dimensionless ionization parameter, $U$, defined as the ratio of the density of ionizing radiation at some position $i$ to the electron density $n_{e}$, is given by

[TABLE]

where $Q$ is the number of ionizing photons per second emitted by the ionizing source, $d_{i}$ is the distance to position $i$ from the ionizing source, and $c$ is the speed of light. In the case of a central AGN, $U_{i}$ decreases as $d_{i}$ increases unless $n_{e}$ also decreases at a faster rate than $1/d^{2}$. The [O III] $\lambda\lambda$4959,5007/[N II] $\lambda\lambda$6548,6584 and [O III] $\lambda\lambda$4959,5007/[S II] $\lambda\lambda$6717,6731 line ratios depend on four parameters: ionization parameter, spectral index (which characterizes the hardness of the ionizing spectrum), gas density, and metallicity. The line ratios will increase outward with an increase in density, a decrease in metallicity, or an increase in ionization parameter.

Since Nuclei B and C have both SDSS and APO/DIS spectra, we can compare their separately measured emission-line ratios to study the ionization gradient. Systematic differences in the absolute flux calibration between SDSS and APO/DIS do not affect line-ratio comparisons. The SDSS spectroscopic aperture is a 3*′′* diameter fiber, whereas the APO/DIS aperture is roughly 1$\farcs$5$\times$2*′′*, i.e., probing a more central region than the SDSS does. As listed in Table 2, the APO 3.5m/DIS measurements of [O III] $\lambda\lambda$4959,5007/[N II] $\lambda\lambda$6548,6584 and [O III] $\lambda\lambda$4959,5007/[S II] $\lambda\lambda$6717,6731 are both larger than those from SDSS. This suggests that the level of ionization decreases with increasing distance to the nucleus for both B and C. The line-ratio gradients cannot be explained by either an outward-increasing density (which is unphysical; the observed electron densities from the SDSS and APO/DIS measurements are consistent within the uncertainties, further ruling this possibility out), or an outward-decreasing metallicity (because the required metallicity would be too high; the line ratios also become insensitive to metallicity at low values, which are expected for both B and C given their low stellar masses). This outward-decreasing ionization parameter is consistent with AGN, in contrast to the outward-increasing ionization parameter expected for distributed ionizing sources such as post-AGB stars, which can produce LINER-like emission-line ratios (e.g., Yan & Blanton, 2012).

We also compare the line ratios [O III]/[N II] and [O III]/[S II] of B against those of A and C (all based on the APO/DIS measurements, which characterize the more central region) to assess whether the gas associated with Nucleus B could be solely ionized by the AGN in Nucleus A and/or C instead of being primarily ionized by its own AGN. The [O III]/[N II] and [O III]/[S II] line ratios of Nucleus B are both comparable to those of Nucleus A and are higher than those of Nucleus C. Based on arguments similar to the above, this comparison suggests that an additional central ionizing source in B is needed and it is not solely ionized by A and/or C. Our conclusion is robust against uncertainties in the dust-reddening correction of the line ratios [O III]/[N II] and [O III]/[S II].

In summary, the line ratios [O III]/[N II] and [O III]/[S II] are most sensitive to the ionization parameter that characterizes the level of ionization. The comparison between the SDSS and APO/DIS-based measurements suggests an outward-decreasing ionization parameter in both B and C, in favor of the AGN hypothesis. The high [O III]/[N II] and [O III]/[S II] line ratios of B as compared against A and C is also in favor of an additional AGN in B, rather than being solely ionized by A and/or C.

3.6 Nature of the Nuclear Ionizing Sources

First, Chandra ACIS-S detects all three nuclei in the soft X-rays (0.5–2 keV) as compact sources whose spatial profiles are consistent with the PSF (Figure 6). The soft X-ray center is consistent with that of the HST $Y$-band nucleus for all three nuclei. Nuclei A and C are both detected as compact (beam size 0$\farcs$19$\times$0$\farcs$19) radio point sources by the VLA in 9.0 GHz (Figure 6). The radio center is consistent with that of the HST $Y$-band nucleus for both nuclei.

Chandra also detects Nucleus A as a compact point source in the hard X-rays (2–8 keV) whose center is consistent with the HST $Y$ band nucleus. Its estimated unabsorbed hard X-ray (2–10 keV) luminosity is 10${}^{42.13^{+0.04}_{-0.05}}$ erg s*-1*, which is comparable to the most X-ray luminous starburst galaxies (Zezas et al., 2001). However, the hard X-ray luminosity cannot be explained by pure star-formation-related processes, because the SFR in the host of A is only modest, estimated as $\sim$8 $M_{\odot}$ yr*-1* (corrected for dust attenuation). The SFR is inferred from HST $U$-band luminosity (corrected for external dust extinction) assuming the empirical calibration based on star-forming galaxies (Hopkins et al., 2003) and verified by an independent estimate based on the 4000 Å break (Brinchmann et al., 2004), which is consistent with the independent constraint ($\lesssim$20 $M_{\odot}$ yr*-1*) based on the VLA radio 9.0 GHz luminosity. On the other hand, the hard X-ray luminosity may be $\sim$10 times higher than our baseline estimate considering systematic uncertainties in the model assumption, making the AGN case even stronger. The high hard X-ray luminosity, particularly when compared to its moderate SFR, unambiguously confirms A as an AGN.

Second, while Nucleus C is undetected in the hard X-rays, it is detected as a compact (beam size 0$\farcs$19$\times$0$\farcs$19) radio source by the VLA in 9.0 GHz (Figure 6). In addition, its SFR is estimated as $\sim$0.1–2 $M_{\odot}$ yr*-1* (corrected for dust attenuation) based on HST $U$-band luminosity and corroborated by the host-galaxy continuum spectral index, which is also consistent with the independent constraint ($\lesssim$2 $M_{\odot}$ yr*-1*) based on the VLA radio 9.0 GHz luminosity assuming that $\lesssim$10% of the radio luminosity comes from a nuclear starburst. This is too low to explain its soft X-ray luminosity by pure star-formation-related processes. Based on the empirical correlation between SFR and LX (Ranalli et al., 2003) observed in star-forming galaxies, the expected soft X-ray luminosity from pure star-formation-related processes is estimated to be L${}^{{\rm SF}}_{0.5-2\,{\rm keV}}{\sim}$1038.7–1040.0 erg s*-1*, which is about an order of magnitude lower than the total soft X-ray luminosity L${}_{0.5-2\,{\rm keV}}{\sim}$10${}^{40.90^{+0.18}_{-0.23}}$ erg s*-1*. Its extinction-corrected $U$-band SFR is more than an order of magnitude lower than that expected from its radio luminosity, assuming it were 100% star formation (i.e., composed of nonthermal synchrotron emission and thermal bremsstrahlung (free-free) emission), according to the empirical correlation between SFR and $L_{R}$ (Bell, 2003) observed in star-forming galaxies.

Finally, while Nucleus B is neither detected in the hard X-rays nor in the radio, its soft X-ray luminosity ($L_{0.5-2\,{\rm keV}}{\sim}$10${}^{41.27^{+0.11}_{-0.13}}$ erg s*-1*) exceeds the luminosity that we would expect from pure star-formation-related processes ($L^{{\rm SF}}_{0.5-2\,{\rm keV}}{\sim}$1039.0–1040.5 erg s*-1*), which is derived from its HST $U$-band and continuum-index inferred SFR estimates ($\sim$0.2–7 $M_{\odot}$ yr*-1* corrected for dust attenuation). Furthermore, the new VLA image allows us to place a stringent radio luminosity upper limit, which independently sets a 3$\sigma$ SFR upper limit of SFR$<$0.8 $M_{\odot}$ yr*-1*. This translates into an even tighter 3$\sigma$ upper limit for $L^{{\rm SF}}_{0.5-2\,{\rm keV}}{<}$1039.6 erg s*-1*, which is more than an order of magnitude lower than the observed total soft X-ray luminosity. The radio-based SFR estimate (and by extension, the soft X-ray luminosity caused by star-formation-related processes) is more robust against uncertainties in the dust-attenuation correction than the $U$-band or host stellar continuum-index estimates. This suggests contribution from an additional excitation source such as an AGN and/or shock heating. In addition, photoionization arguments suggest that its Seyfert-type BPT diagnostic line ratio is unlikely to be caused by pure star-formation-related process and/or shock heating, or primarily by the AGN in Nuclei A and/or C. The line ratios [O III]/[N II] and [O III]/[S II] are most sensitive to the ionization parameter that characterizes the level of ionization. Comparison between the SDSS (from a 3*′′* diameter aperture) and APO/DIS-based (roughly from a 1$\farcs$5$\times$2*′′* aperture, i.e., probing a more central region than the SDSS) line ratios suggests an outward-decreasing ionization parameter in both B and C, in favor of the AGN hypothesis. The high [O III]/[N II] and [O III]/[S II] line ratios of B as compared against A and C are also in favor of an additional AGN in B, rather than being solely ionized by A and/or C. Finally, while shocks could also produce excess of soft X-ray emission, the Seyfert-type optical emission-line ratio is difficult to accommodate in a purely shock-heating scenario.

In summary, our complementary multiwavelength follow-up observations lend further support to the triple BH scenario suggested by the optical slit spectroscopy. Taken together, our comprehensive observations strongly suggest that all three nuclei are AGN, making SDSS J0849+1114 the most unambiguous case known to host a kpc-scale trio of MBHs.

4 Discussion

4.1 Black Hole Mass

We estimate the BH mass for each nucleus inferred from the host-galaxy total stellar mass and bulge stellar velocity dispersion assuming the $M_{\bullet}$–$M_{\ast}$ and $M_{\bullet}$–$\sigma_{\ast}$ relations (Kormendy & Ho, 2013). For the $M_{\bullet}$–$M_{\ast}$ relation, we prefer the host-galaxy total stellar mass over the bulge stellar mass because (i) there is empirical evidence that the correlation between BH mass and total stellar mass is tighter than that between BH mass and bulge stellar mass in low-redshift active galaxies (Kormendy & Ho, 2013), and (ii) our bulge and disk decomposition may be uncertain, making the total stellar mass estimate more robust than the bulge stellar mass estimate. We adopt the relation given by

[TABLE]

which is calibrated based on a sample of 262 local broad-line AGN with virial BH mass estimates (Reines & Volonteri, 2015). The rms scatter in the relation is 0.55 dex, which includes the measurement error in the virial BH mass (0.50 dex) and the best-fit intrinsic scatter (0.24 dex). The stellar masses in the calibration sample of Equation 7 are based on the M/L ratios in Zibetti et al. (2009), which assumes a Chabrier IMF. Our stellar mass uses the M/L ratios in Bell et al. (2003) assuming the scaled Salpeter IMF. The Bell et al. (2003) masses are systematically higher than the Zibetti et al. (2009) masses by 0.21 dex (Reines & Volonteri, 2015). To correct for this systematic offset, we subtract 0.21 dex from our stellar mass estimate as the input for Equation 7. Table 1 lists the estimated BH mass for each nucleus from the $M_{\bullet}$–$M_{\ast}$ relation.

For the $M_{\bullet}$–$\sigma_{\ast}$ relation, we adopt the relation given by

[TABLE]

which is calibrated based on a sample of 71 local active galaxies with virial BH mass estimates (Greene & Ho, 2006). The intrinsic scatter in the relation is $<$0.4 dex. Table 1 lists the estimated BH mass for each nucleus from the $M_{\bullet}$–$\sigma_{\ast}$ relation, which is consistent with the estimate from the $M_{\bullet}$–$M_{\ast}$ within the uncertainties, lending support to our BH mass estimates. While still broadly consistent with the stellar-mass-based estimates, the stellar velocity-dispersion-based estimates may result in overestimated BH masses, considering that stellar velocity dispersions are systematically higher in mergers than in isolated AGN host galaxies (e.g., Liu et al., 2012). Our estimates are also broadly consistent with independent constraints assuming the fundamental place relation of BH accretion (Merloni et al., 2003; Gültekin et al., 2019).

4.2 Merging Timescales

We use Chandrasekhar dynamical friction timescale arguments to estimate characteristic timescales for the galaxy-galaxy mergers and the inspiraling timescales of the BHs. Assuming singular isothermal spheres for the host-galaxy density profiles and circular orbits, we estimate the dynamical friction timescale for a satellite galaxy with velocity dispersion $\sigma_{{\rm S}}$ inspiraling from radius $r$ in a host galaxy with velocity dispersion $\sigma_{{\rm H}}$ as (Chandrasekhar, 1943)

[TABLE]

where $\Lambda$ is the Coulomb logarithm with typical values of 3${\lesssim}\ln\Lambda{\lesssim}$30 (Binney & Tremaine, 1987). $\ln\Lambda{\sim}$2 for equal-mass mergers (Dubinski et al., 1999). We assume $\ln\Lambda{\sim}$6. With A as the host and B and C as satellites, B and C will merge with A in $t_{{\rm fric}}{\sim}$56${}^{+36}_{-28}$ Myr and $\sim$49${}^{+77}_{-25}$ Myr, corrected for projection effects assuming random orientation.

After merging together with their individual hosts, the three BHs will in-spiral to the center of the merged galaxy under dynamical friction with the stellar background. Assuming a singular isothermal sphere for the density distribution of the merged host galaxy, the dynamical friction timescale of a BH of mass $M_{\bullet}$ on a circular orbit of radius $r$ is estimated as

[TABLE]

where $\ln\Lambda{\sim}$6 for typical values (Binney & Tremaine, 1987). We estimate $t_{{\rm fric}}{\sim}$1.9${}^{+2.0}_{-0.9}$ Gyr and 1.5${}^{+0.9}_{-1.0}$ Gyr for the BH in B and C to reach the center of the merger products to form a gravitationally bound triple BH system with A, assuming a radius of $r$=1 kpc at the start of the inspiraling.

We have neglected the effects of friction by gas, which will accelerate the merger processes. On the other hand, we have neglected the effects of tidal stripping of stars in the satellites, which will delay the merger. We have also neglected the effect of BH accretion, which may change the masses of the BHs and the resulting merger timescales. The actual merger involving three stellar and BH components with gas friction and gas accretion onto three BHs is likely much more complicated. Numerical simulations with model parameters tailored to those of SDSS J0849+1114 are needed to make more realistic predictions of the subsequent merger of the triple system (Renaud et al., 2010). Regardless of the actual sequence of the mergers, the BHs are likely to form a gravitationally interacting triple system if the orbit of the third MBH decays rapidly enough before the first two BHs coalesce.

5 Summary and Future Work

We have reported the discovery of SDSS J0849+1114, the first case of a triple Type 2 Seyfert nucleus. It represents three active BHs hosted in a three-way galaxy merger. The target was selected from the largest sample of optically selected candidate AGN pairs based on the SDSS data. However, the archival SDSS data only have spectra for two of the three nuclei; the SDSS fibers were also too large to clearly separate the emission from the nucleus. To unambiguously determine the nature of the triple nucleus, we have conducted a comprehensive multiwavelength follow-up campaign. We summarize our main findings as follows:

• •

First, by conducting new spatially resolved optical spectroscopy, we have classified all three nuclei as Type 2 Seyferts based on the classical BPT diagram. This establishes J0849+1114 as the first known case of a triple Type 2 Seyfert nucleus.

• •

Furthermore, to confirm the excitation mechanism of the triple Type 2 Seyfert nucleus, we have further conducted a complementary follow-up campaign including Chandra X-ray imaging spectroscopy, HST $U$- and $Y$-band imaging, and VLA radio imaging. These new comprehensive multiwavelength observations strongly suggest that J0849+1114 hosts a kpc-scale triple AGN. This represents the most robust evidence for three active MBHs in the process of merging. The discovery provides a unique verification of the hierarchical MBH-formation paradigm with important implications for the detection of low-frequency gravitational waves.

• •

By modeling the host-galaxy photometry and internal dust extinction of SDSS J0849+1114, we have estimated the stellar masses of the three merging components to be $\sim$10${}^{11.3}M_{\odot}$, 10${}^{10.0}M_{\odot}$, and 10${}^{10.5}M_{\odot}$ for A, B, and C. Assuming the empirical correlation between BH mass and host total stellar mass observed in local broad-line AGN (Reines & Volonteri, 2015), the inferred BH masses are $\sim$10${}^{7.5}M_{\odot}$, 10${}^{6.4}M_{\odot}$, and 10${}^{6.7}M_{\odot}$, consistent with independent estimates based on the host-galaxy stellar velocity dispersion within the uncertainties.

SDSS J0849+1114 provides a unique verification of the hierarchical merger-tree model for MBH formation (e.g., Volonteri et al., 2003). Smaller galaxies merge to form larger ones (Toomre & Toomre, 1972). After galaxies merge, massive binary BHs should form through dynamical friction (Begelman et al., 1980; Milosavljević & Merritt, 2001; Yu, 2002). When a third galaxy enters before the first two BHs coalesce, an MBH triple should form (Valtonen, 1996). Triple MBHs are important for understanding galaxy formation. They could scour out stellar cores much larger than those formed around binaries (Hoffman & Loeb, 2007). In addition, triple MBHs represent a unique laboratory for studying the chaotic dynamics of general relativity three-body interactions (Blaes et al., 2002; Merritt, 2006). Their end products may include a merger of all three BHs, formation of a highly eccentric binary, or even ejection of three free BHs (Lousto & Zlochower, 2008). The inner binary may have a very high eccentricity, releasing intense bursts of gravitational waves sought by pulsar-timing arrays and the Laser Interferometer Space Antenna (Amaro-Seoane et al., 2010). Three-body interactions also represent one of the main mechanisms to solve the final-parsec problem for rapid binary coalescence (Valtonen, 1996; Blaes et al., 2002; Hoffman & Loeb, 2007).

Despite significant merit and intense interest, direct evidence for triple MBHs is still lacking. Their separations (less than a few parsecs) are too small to resolve beyond the local universe with current facilities. Kpc-scale triple AGN represent a promising precursor of compact triple MBHs. Using a simple stellar dynamical friction argument (§4.2), we have estimated that the trio in SDSS J0849+1114 may form a bound MBH triple in $\lesssim$2 Gyr.

Future much deeper X-ray observations with Chandra may detect and resolve the two fainter nuclei in the hard X-rays and characterize their X-ray spectral properties. Deeper VLA observations, now under way, may detect the faintest Nucleus B in the radio and measure its radio spectral index. SOFIA/FORCAST imaging may measure the mid-IR spectral energy distributions and marginally resolve the three nuclei. Finally, high spatial resolution, integral field unit spectroscopy in the optical and/or near-IR may map out the ionization and velocity fields around the triple AGN.

The selection of SDSS J0849+1114 is limited because its parent galaxy sample is based on the SDSS. Future large near-IR spectroscopic galaxy surveys (Takada et al., 2014; The MSE Science Team et al., 2019) may uncover a larger population of galaxies harboring triple AGN in the higher redshift universe where mergers are believed to be more frequent.

Acknowledgements

We thank S. Gilmore, M. Claussen, and N. Ferraro for help with our VLA observations, and the anonymous referee for a prompt and helpful report. X.L. acknowledges support by NASA through Chandra Award Number GO3-14104X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS 8-03060. Support for Program number HST-GO-13112 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. M.H. and Z.L. acknowledge support by the National Key Research and Development Program of China (2017YFA0402703). M.K. acknowledges support by NSFC Youth Foundation (No. 11303008) and by the Astronomical Union Foundation under grant No. U1831126. Y.S. acknowledges partial support from an Alfred P. Sloan Research Fellowship and NSF grant AST-1715579. This work makes extensive use of SDSS-I/II and SDSS-III/IV data (http://www.sdss.org/ and http://www.sdss3.org/). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This paper includes data gathered with the Apache Point Observatory 3.5 m telescope, which is owned and operated by the Astrophysical Research Consortium.

This research has made use of software provided by the Chandra X-ray Center in the application packages CIAO, ChIPS, and Sherpa.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group,Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

Facilities: APO 3.5m (DIS), CXO (ACIS), HST (WFC3), Sloan, VLA.