TXS 2116−077: A Gamma-Ray Emitting Relativistic Jet Hosted in a Galaxy Merger

We obtained a high-resolution J-band image of TXS 2116−077 with the InfraRed Camera and Spectrograph (IRCS; Kobayashi et al. 2000), in combination with its AO system AO 188 (Hayano et al. 2010), mounted at the Subaru Telescope on 2018 June 16. We used a resolution of 0052 per pixel for a corresponding 54'' field of view. The Laser Guide Star mode was used for AO and observations were performed in imaging mode. The total integration time was 34 minutes while accumulating 17 frames with an exposure time of 120 s per frame using a 5 point dither pattern. The natural seeing was ∼07 and the AO-corrected seeing was 02.

The IRCS data was analyzed with Image Reduction and Analysis Facility (IRAF; Tody 1986). This includes the creation of a bad pixel mask, background subtraction, flat-fielding estimate of the dither offsets, and final stacking of the data. A nearby star, SDSS J211853.33−073214.3 (R.A.: 21^h18^m5333, decl.: −07°32'1434), located ∼14'' northeast of the γ-NLSy1 was used as a reference for the flux calibration.

ISIS at WHT. Long-slit spectrophotometric observations of TXS 2116−077 were obtained in service time with the Intermediate-dispersion Spectrograph and Imaging System (ISIS) double-arm spectrograph mounted on WHT of the Isaac Newton Group at the Roque de los Muchachos Observatory on the Spanish island of La Palma. They were acquired on 2018 August 13 with an average seeing of 09. To cover both nuclei detected in the Subaru image, we used a slit position angle of 1167. Grating R158B was used in the blue arm centered at 4500 Å with a wavelength range between 3370 and 5400 Å and 1.45 Å pixel⁻¹ dispersion. In the red arm, grating R158R centered in 7500 Å covered between 5300 and 10000 Å with a dispersion of 1.82 Å pixel⁻¹. The slit width was 12. A total of six exposures of 1800 s each were taken simultaneously in each arm. The reduction was performed using IRAF routines in the usual manner. For both arms, the standard star G191-B2B was used to produce final sky subtracted and flux calibrated spectra.

OSIRIS at GTC. Long-slit spectrophotometric observations were obtained under the Director's Discretionary Time (DDT) program using the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS) spectrograph mounted at GTC located at the Roque de los Muchachos Observatory on the island of La Palma, Spain. They were taken on 2018 August 16 with an average seeing of 09 during dark time using grating R2000B and a slit width of 08. The wavelength range covered between 3950 and 5650 Å at a dispersion of 0.86 Å pixel⁻¹. A total of four exposures of 1800 s each were taken. The reduction was performed using IRAF routines following the standard procedures. The standard star Ross 640 was used to obtain the sensitivity curve and perform the flux calibration.

Multi-Espectrógrafo en GTC de Alta Resolución para Astronomía (MEGARA) is an optical IFU and multi-object spectrograph installed at the GTC (Gil de Paz et al. 2016). The MEGARA IFU encompasses 567 contiguous hexagonal fibers, each with a long diagonal of 062, resulting in a 12.5 × 11.3 arcsec² field of view in the shape of a rectangle. Sky observations are obtained simultaneously with eight fiber bundles with seven dedicated hexagonal fibers each located at the outermost parts of the field of view (R > 15). MEGARA IFU spectra of TXS 2116−077 were obtained under DDT program on the night of 2018 October 1 using the LR-I VPH grating that covers from ∼7200 to 8650 Å with a reciprocal linear dispersion of 0.37 Å pixel⁻¹. We obtained six exposures of 1690 s each. The seeing conditions were around 06. To flux-calibrate the LR-I spectra, several exposures of the flux standard star HR 7950 were obtained under similar airmass and seeing conditions. We also obtained calibration data, including ThNe and ThAr arc lamps and halogen lamps.

The raw IFU spectra were processed using the MEGARA data reduction pipeline¹² (MDRP; S. Pascual 2020, in preparation). The main steps include: bias correction, trimming, identification of the position of the spectra on the detector along the dispersion axis, model mapping and extraction of each individual spectrum, fiber-flat correction, wavelength calibration, flux calibration, and sky subtraction. The final products of the MDRP are row-stacked spectra (RSS) 2D images containing 623 fiber spectra. The Hα line shows a second, broad component in some spaxels around the Seyfert 1 nucleus. Thus, we used a multi-component Gaussian model to fit the [N ii] + Hα window. For the sake of visualization, we have transformed the hexagonal spaxel grid of the RSS data into a square pixel shape grid of 025 per pixel.

The archival Very Large Array (VLA) data of TXS 2116−077 at 8.4 GHz, observed on 1998 May 18, was analyzed using the Astronomical Image Processing System (AIPS), following standard procedures. The flux scale was set using the standard flux calibrator 3C 286. The flux and phase calibrations were transferred to the target via the secondary calibrator. A few rounds of self-calibration were performed to correct for phase variations, which improved the image quality. The resolution of the final image is 032 × 023 at a position angle of −14°. The root mean square (rms) noise floor of the map is 42 μJy beam⁻¹. The γ-NLSy1 remained unresolved at this frequency and resolution. We used the Gaussian fitting program JMFIT in AIPS to get the flux density. The peak and integrated flux density is nearly identical confirming the compact nature of the source. The integrated flux density of the target is 91.87 mJy. To compute the error on the flux, we have added the calibration error, Gaussian fitting error, and the rms noise of the map in quadrature, giving the final error as 0.17 mJy. The companion Seyfert 2 galaxy, on the other hand, remain undetected and we derived the 3σ flux upper limit as 0.13 mJy.

A deep radio imaging of TXS 2116−077 with Very Large Baseline Array (VLBA) has recently been performed on 2018 March 30. We analyzed this publicly available VLBA observation at 6 cm (5 GHz). A sustained data-recording rate of 512 Mbits s⁻¹ in two-bit sampling was used. The frequency band was split into four intermediate frequencies (IFs) of 64 MHz bandwidth each, for a total synthesized bandwidth of 256 MHz. Each IF was in turn split into 256 channels of 250 KHz bandwidth. The data were corrected for Earth-orientation and ionospheric effects. The source BL Lac was used as bandpass calibrator. The data were correlated at the NRAO data processor using an averaging time of 2 s. We performed standard a-priori gain calibration using the measured gains and system temperatures of each antenna. This calibration, as well as the data inspection and editing, were done within AIPS. The flux scale was set using the calibrator BL Lac. The observation resulted in a beam of 5.3 × 2.4 millarcsec², at a PA of −2°. A compact radio emission positionally consistent with the γ-NLSy1 galaxy is found with no hints of any extended features or knots. The peak flux maximum is located at the position the R.A.: 21^h18^m5296 and decl., J2000 = −07°32'2758, with a flux density of 27.93 ± 0.23 mJy b⁻¹. There is no detection of the companion Seyfert 2 galaxy, with an upper limit of 0.13 mJy beam⁻¹.

Chandra X-ray observatory carried out a DDT observation of TXS 2116−077 on 2018 September 15 for a net exposure of 10 ksec. We used the nominal aim point of the Advanced CCD Imaging Spectrometer chip (ACIS-S3) and the data were collected in a full readout mode resulting in 3.241 s readout time. The data were reprocessed using chandra_repro and an exposure-corrected image was generated using fluximage. We scanned the Chandra image with the source detection algorithm wavedetect to search for X-ray emitting objects. A bright X-ray source (R.A., J2000 = 21^h18^m5295 and decl., J2000 = −07°32'2754) coincident with the γ-NLSy1 is identified with the separation between the X-ray and radio positions as 019. The excellent spatial resolution of Chandra (∼05) allowed us to resolve the merging system, however the companion Seyfert 2 galaxy remains below the detection limit.

We used the task specextract to extract the source and background spectra and corresponding ancillary and response matrix files. A source region of 15 centered at the γ-NLSy1 galaxy was adopted and we selected the background as a circle of 10'' from a nearby source-free region. We rebinned the source spectrum to have 1 count per bin and performed the fitting in the energy range of 0.5−7 keV using C-statistics (Cash 1979) in XSPEC (Arnaud 1996).

The results acquired from the Subaru imaging and subsequent follow-up observations of TXS 2116−077 confirm its ongoing merger with a Seyfert 2 galaxy. This is probably the first unambiguous detection of a γ-ray emitting jet hosted by a merging system and hints that relativistic jets can be triggered by galaxy mergers (Chiaberge et al. 2015). Note that a few recent studies have also reported the identification of minor merging events associated with NLSy1 galaxies (e.g., Järvelä et al. 2018; Olguín-Iglesias et al. 2020), however, due to a lack of spectroscopic redshift measurement of the companion, strong claims could not be made. Berton et al. (2019) presented the spectroscopic confirmation of the merger of a radio-loud NLSy1 galaxy IRAS 20181−2244 (z = 0.185), however, this object is yet to be detected in the γ-ray band and thus its extent of relativistic beaming is uncertain. Furthermore, the small physical separation (∼12 kpc) and the presence of AGN activity in both galaxies indicates that this is a merger in the late stage of evolution, i.e., close to coalescence (Di Matteo et al. 2005; Van Wassenhove et al. 2012). Both galaxies are embedded in a common, 70 kpc extended envelope composed of gas and stars, and the companion galaxy exhibits morphological features associated with tidal disruption (Figure 1).

The photometric decomposition of the Subaru image was performed with the 2D fitting algorithm GALFIT (Peng et al. 2010). The star SDSS J211853.33−073214.3 was used to model the point-spread function (PSF) and remove the AGN contamination. A nearby source-free region was used as the sky background. First, we fitted both galaxies only with the PSF and then added components, mainly Sérsic profiles, based on the leftover emission in the residual image and the reduced χ² value. Typically, a smaller value of the Sérsic index (n ≲ 2) is associated with galaxies with disk-like morphology and pseudobulges, and larger values of n(>2) with elliptical galaxies and classical bulges.

The best-fitted model includes two Sérsic functions for each AGN and three Fourier components (see Peng et al. 2010, for details) to take into account the extended emission. The resultant fitting parameters are provided in Table 1 and the results of the fitting are shown in Figure 5. For both γ-NLSy1 and the companion, first Sérsic functions with n = 1.48 ± 0.06 and 0.73 ± 0.03, respectively, resemble that of pseudobulges in late-type galaxies. The physical interpretation of the second Sérsic components is tedious due to disturbed morphologies of both galaxies. However, the small Sérsic index for this component indicates the presence of a disk/bar in both galaxies. Altogether, it can be concluded that both merging galaxies are "young" and host pseudobulges, similar to what is commonly known for other NLSy1 (Mathur 2000; Järvelä et al. 2018; Olguín-Iglesias et al. 2020). This reinforces the hypothesis that γ-NLSy1 galaxies, and NLSy1s in general, are young sources within the AGN evolution scheme.

Zoom In Zoom Out Reset image size

The BPT diagram (Baldwin et al. 1981) is a powerful tool to explore the excitation mechanism of the ionized gas. We used the WHT and GTC long-slit spectroscopic observations and extracted spectra at eight locations along the slit with the motivation to understand the physical behavior of the merging environment. In Figure 6, we show the location of seven of the eight spectra extracted¹³ from the long-slit spectroscopy (large red dots) over the background of the BPT diagram of Kewley et al. (2006) for 85,224 galaxies from the SDSS data release 4. Because these seven spectra are extracted from the long-slit along a stretch of 15'' (see black lines in Figure 3), we can firmly state that a large fraction of the gas along the slit in the merging system is affected by AGN and shocks, i.e., physical processes other than due to hot, massive stars. The BPT diagram was also used to elucidate the nature of the enhanced Hα emitting ring surrounding the γ-NLSy1 galaxy observed in the MEGARA IFU image (Figure 3, top panel) which confirmed that it is excited by the merger induced shocks and the AGN emission.

Zoom In Zoom Out Reset image size

After correcting for 3.6 Å instrumental resolution, the full width at half maximum (velocity dispersion) value for the emission line [O iii] 5007 Å is 10.29 Å (208 km s⁻¹) for the Seyfert 2 nucleus. The FWHM[O iii] can be converted to the stellar velocity dispersion (σ) by σ = FWHM[O iii]/2.35 (Nelson & Whittle 1995; Wang & Lu 2001). Using the known statistical correlation between the mass of the central black hole and σ (Wang & Lu 2001; Grupe & Mathur 2004; Gültekin et al. 2009), this value transforms to a black hole mass of ∼4 × 10⁶ M_⊙ for the companion Seyfert 2 galaxy. For the γ-NLSy1 galaxy, on the other hand, we used the well calibrated empirical relations (see Shen et al. 2011) adopting the conventional Hβ FWHM (∼1900 km s⁻¹) and continuum luminosity at 5100 Angstrom (∼10⁴⁴ erg s⁻¹; see also Rakshit et al. 2017). The derived black hole mass is ∼3 × 10⁷ M⊙, similar to that reported in literature (e.g., Rakshit et al. 2017). Note that the typical uncertainties associated with these approaches are 0.3−0.5 dex (Vestergaard & Peterson 2006; Gültekin et al. 2009; Shen et al. 2011). Furthermore, since the black hole mass and the host galaxy properties are significantly correlated (Di Matteo et al. 2005; Reines & Volonteri 2015), the total stellar mass is estimated as ∼1 × 10¹¹ M_⊙ and ∼2 × 10¹⁰ M_⊙ for the γ-NLSy1 and the companion galaxy, respectively. Based on the large ratio of the stellar masses, the observed event can be classified as a "minor" merger (Lotz et al. 2011).

We use the spectral population synthesis code starlight (Cid Fernandes et al. 2005) to derive key population parameters, such as mean age and metallicity of the merging system. The continuum stellar spectra (masking the emission lines) were fitted for each of the eight spectra extracted from the spatial locations along the OSIRIS and ISIS slits, as demonstrated in the left panel of Figure 7. We retrieved the star formation history, the dust extinction, mass weighted stellar metallicity, and luminosity weighted stellar ages and show them in the right panel of Figure 7. We used a set of 163 single stellar populations sampling the 3 Myr to 9.75 Gyr¹⁴ age range with 33 logarithmically spaced age bins and five metallicities between 0.02 and 2.5 solar (S. Charlot and G. Bruzual 2007, private communication), assuming a Chabrier initial mass function. Further details about starlight can be found on the project web page http://www.starlight.ufsc.br, and a review of stellar populations fitting techniques is available at http://www.sedfitting.org (Walcher et al. 2011).

Zoom In Zoom Out Reset image size

As can be seen in Figure 7, the results from the stellar population synthesis calculations performed on the GTC and WHT optical spectra show average (luminosity weighted) stellar ages in the range ∼0.5–2.5 Gyr, older in the outer halo and younger in the inner region around the Seyfert 2 nucleus. The presence of the younger stellar population in proximity to the merging nuclei indicates that the starburst activity has recently begun in the central regions. This result is expected if the spectra are dominated by merger driven activity in the past ∼0.5–2 Gyr (Moreno et al. 2015). This is because the star formation history that develops along the evolution of a merging pair of galaxies depends on the distance between the galaxies during the different binding passages. For example, Figures 3 and 4 in Moreno et al. (2015) show the star formation evolution from the approach and first passage until the coalescence. Clearly, the window of higher star formation, i.e., when more gas is available to form stars and also to fall into the nuclei, is in the approximate age range 0.5−2 Gyr. The detailed history of star formation and nuclear activity also depends on the relative masses and morphology of the two merging galaxies and on the relative angular momenta of the spirals and the merging orbit; but models show that nuclear activity follows after the main bouts of star formation. Indeed, Van Wassenhove et al. (2012) find that strong dual AGN activity occurs during the late phases of the mergers, at small separations (≲10 kpc). Thus, the results from the analysis of the stellar populations and the nuclear activity in TXS 2116−077, especially the merger timescale, match quite precisely the expectations from the numerical simulations.

From the total narrow Hα flux measured in the MEGARA data (15.7 mJy), the star formation rate in the system is 3.8 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$, using the established correlation of star formation rate with Hα emission (Kennicutt 1998). Using the correlation between the global star formation rate and the galaxy stellar mass (González Delgado et al. 2016; Cortijo-Ferrero et al. 2017), this places TXS 2116−077 just within but somewhat high in the star-forming main sequence, implying that the star formation is possibly enhanced by the merging process when compared to normal late-type spiral galaxies (González Delgado et al. 2016).

Since the radio jet in TXS 2116−077 remains compact down to milliarcsecond scale, we have computed the upper limit to its projected length at z = 0.26 as 0.02 kpc, using the VLBA beam size of (5.3 × 2.4) mas². Considering a viewing angle of 3°, typical for γ-ray emitting beamed AGNs (Ghisellini et al. 2014; Paliya et al. 2017), the true length turns out to be <1 kpc after correcting for projection effects. This gives the kinematic age of the jet as <15 kyr, assuming a typical jet velocity of 0.1 c (e.g., Gitti et al. 2013). The age will be even shorter for the larger jet velocity/viewing angle. This finding suggests a "young" jet in TXS 2116−077, which is in line with the hypothesis that NLSy1 jets are likely to be in the early phase of evolution (Singh & Chand 2018). Note that the radio emission in NLSy1 galaxies could be originated from the star formation activities and/or due to nonrelativistic kiloparsec-scale outflows (Gallimore et al. 2006). However, the observation of the flux variability, especially in the γ-ray band, and also the detection of a flat GHz frequency spectrum (Yang et al. 2018) indicates that the radio emission detected from TXS 2116−077 is jet dominated and hence the derived kinematic age is reliable.

A comparison of the jet kinematic age (<15 kyr) with the merger timescale (0.5−2 Gyr) suggests the jet to be considerably younger than the merger. This observation provides supportive evidence for a scenario in which relativistic jets are triggered by galaxy mergers. Though high-quality observations, similar to those presented here, are not available for many jetted AGNs, recent observing campaigns focused on the high-resolution host galaxy imaging of radio-loud NLSy1 galaxies are also converging to this scenario (see, e.g., Berton et al. 2019; Olguín-Iglesias et al. 2020). In particular, Olguín-Iglesias et al. (2020) studied a sample of γ-NLSy1s and reported that minor mergers could be responsible for the observed nuclear activity in this class of AGN. Though host galaxy observations of a few other sources revealed a rather secular driven growth of the system (Kotilainen et al. 2016) and a few studies also claimed elliptical hosts (e.g., D'Ammando et al. 2017), the majority of the γ-NLSy1s undoubtedly reside in young, late-type galaxies in which jet launching is likely to be triggered via galaxy mergers. Similar results are found in various studies focusing on more powerful, radio-loud quasars (Ramos Almeida et al. 2013; Chiaberge et al. 2015; Olguín-Iglesias et al. 2016). The fact that radio-loud NLSy1s are likely to be living in denser large-scale environments (Järvelä et al. 2017) where mergers could be frequent, also provides supportive evidence for the merger triggered jet launching. Combining this with the idea of the young nature of NLSy1s, it can be concluded that γ-NLSy1 galaxies represent the beginning phase of the jet activity.

There is significant support of the idea that rotating supermassive black holes produce jets and their origin is connected to magnetic fields generated in accretion disks of spinning black holes (Blandford & Znajek 1977; Meier 1999). Recent 3D relativistic magnetohydrodynamic simulations support this spin paradigm (see, e.g., Tchekhovskoy et al. 2011) and demonstrate that prograde disk accretion onto a rapidly spinning black hole most efficiently generates jets. Combining this with the idea that galaxy mergers trigger AGNs (Di Matteo et al. 2005; Hopkins & Quataert 2011) leads to the expectation that a merger may produce the needed high spin. Galaxy mergers can efficiently drive gas inflows all the way to the central regions of the galaxy and all the way to the accretion disk that feeds the central black hole.