© B. Mazzilli Ciraulo et al. 2021

1. Introduction

The importance of galaxy interactions and merger events is nowadays beyond dispute as hierarchical galaxy evolution is widely adopted and multiple observational as well as theoretical works have been carried out to investigate the various mechanisms related to the evolution of interacting galaxies (e.g. van Dokkum et al. 2015; Cappellari 2016). The improvement in data resolution has made it possible to highlight various observational signatures of mergers, such as colour change (Alonso et al. 2012), morphological disruption (Ellison et al. 2010; Casteels et al. 2013), star formation enhancement (Patton et al. 2013), and merger-induced nuclear activity at low redshift for optical and mid-IR active galactic nuclei (AGNs; Ellison et al. 2019). The presence of double nuclei (such as in NGC 3526; Sakamoto et al. 2014) is a clear signature of mergers as the existence of tidal tails (Mesa et al. 2014) and, alternatively, spectral features, such as double-peaked profiles in molecular gas emission lines, have also been interpreted as merger signatures (e.g. Greve et al. 2005; Weiß et al. 2005). While the timescales of the merging are well monitored statistically through numerical simulations (e.g. Somerville & Davé 2015; Nevin et al. 2019) and for individual objects (e.g. Lahén et al. 2018), dating the epochs of these observed processes remains a challenge. Lastly, some important questions remain unanswered regarding the possible impact of AGNs on the merging stage (e.g. Alonso et al. 2007; McAlpine et al. 2020; Shah et al. 2020; Sharma et al. 2021) and whether the quenching is related to the merging (e.g. Ellison et al. 2018; Hani et al. 2020; Das et al. 2021).

In this article we present a multi-wavelength analysis of a galaxy merger (J221024.49+114247.0) at z = 0.09, identified through the detection of double-peaked emission-line profiles in the central spectrum of the Sloan Digital Sky Survey (SDSS; Strauss et al. 2002). Two gas and stellar counter-rotating components have been identified and analysed. In Sect. 2 we discuss how we extract properties from the Mapping Nearby Galaxies at Apache Point Observatory (MaNGA; Bundy et al. 2015) data cube to identify the two components. In Sect. 3 we discuss the molecular gas data, as well as other archives (Very Large Array) and properties used from value-added catalogues. In Sect. 4 we analyse the properties of this merger. In Sect. 5 we discuss our results. We adopt throughout this work a Λ cold dark matter (CDM) cosmology with parameters Ω_Λ = 0.7, Ω_M = 0.3, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. MaNGA data analysis

The galaxy J221024.49+114247.0 is one optical merger identified in a sub-sample (Mazzilli Ciraulo et al., in prep.) obtained by cross-matching the summary file from the MaNGA Data Reduction Pipeline (DRP; SDSS Data Release DR14, Abolfathi et al. 2018) and a double-peak galaxy catalogue produced by Maschmann et al. (2020), relying on the Reference Catalogue of Spectral Energy Distributions (RCSED) from Chilingarian et al. (2017). The basic properties of this source, provided by different value-added catalogues, are given in Table 1. The Legacy survey (Dey et al. 2019) reveals a disrupted object, and the integrated spectrum within the 3 arcsec SDSS fibre shows double-peaked emission lines. The effective point spread function (PSF) of the MaNGA data is about 2.5 arcsec and corresponds to a spatial resolution of 4.3 kpc at the redshift of the source (z = 0.09228), and the data are gridded on 0.5 × 0.5″ spaxels.

We produced 2D maps based on the derived properties from the MaNGA Data Analysis Pipeline (DAP; SDSS-IV, Westfall et al. 2019). The main emission lines peak in a compact central region. The velocity field does not display any standard rotating motion and appears to be strongly irregular, with values covering a very wide range (−525 km s⁻¹ < V_gas < 445 km s⁻¹). As displayed in Fig. 1, the velocity dispersion is also remarkable, with spaxels exhibiting values greater than 400 km s⁻¹ and not centred. The double-peak features of these spaxels are not well fitted by the single component fit, and some of them are masked during the data processing. In order to further investigate this complex system, we developed an optimised multi-component analysis.

Fig. 1. Velocity dispersion of the Hβ emission line derived by the MaNGA DAP. The MaNGA field of view is represented as a grey footprint. The included panel shows a portion of the continuum-subtracted spectrum from one spaxel where the model (in orange) does not manage to fit the spectral features.

2.1. Gas kinematics derived from the multi-component approach

We developed a fitting procedure, adjusting two components for the selected emission lines. We ran and applied it to the Voronoi binned continuum-subtracted spectra and fitted Hβ, [OIII] doublet, Hα, and [NII] doublet lines. As in the MaNGA DAP, we assigned the same velocity and velocity dispersion to all these lines. We were able to highlight the presence of two components along the line of sight. We first selected reliable spaxels, relying on a signal-to-noise ratio (S/N) threshold of 3. When the S/N criterion is not met by the two fitted components, we kept the single-Gaussian results and associated them with the main galaxy. For the secondary component, we based our spaxel selection on an F-test (Mendenhall & Sincich 2011) to ensure that the double-Gaussian fit was better than the single Gaussian one, that a velocity difference between both peaks was significantly greater than the MaNGA data spectral resolution, and that there was an amplitude ratio between 1/3 and 3 between both peaks to guarantee that one of them was not completely suppressed.

We display the derived properties for each of them on the 2D maps shown in Figs. 2 and 3. Figure 2 displays the extinction-corrected Hα line fluxes, the extinction derived from the Balmer decrement, and the gas metallicity using the O₃N₂ calibrator, while Fig. 3 exhibits the gas kinematics of the two components with their velocity and velocity dispersion. The spatially resolved maps confirm the contribution of two separate components: The first is related to a galaxy with a fairly regular velocity field, whereas the second is detected in a smaller region of the field of view, shows a slightly stronger Hα line flux, and shows a velocity field counter-rotating with respect to the main object.

Fig. 2. Gas properties derived from our multi-component approach. First and second rows: maps for the first and second components detected in MaNGA data, respectively. First column: Hα extinction-corrected flux (in erg s−1 Å−1 cm−2 per spaxel), second column: extinction computed from the Balmer decrement, and last column: oxygen gas-phase abundance derived using the O3N2 calibrator. The black crosses indicate the position of the extinction-corrected Hα flux peak for the represented component. The MaNGA PSF is displayed as a hatched grey circle in the bottom-left corner of the panels.

Fig. 3. Gas and star kinematics. First row: maps for the first component, and second row: second component detected in MaNGA data. First and second columns: respective ionised gas and stellar velocity fields (in km s−1). The dotted black lines refer to the computed position angles. Third and fourth columns: respective velocity dispersion fields (in km s−1) for the gas and the stars. As in Fig. 2, the black crosses indicate the position of the extinction-corrected Hα flux peak for the represented component.

2.2. Multi-component analysis of the stellar kinematics

To explore whether the double-peaked structure in the emission lines has associated counterparts in the stellar kinematics, we applied an analysis workflow used for the study of galaxies that host counter-rotating stellar discs (Katkov et al. 2013, 2016). Firstly, we fitted spectra using a full spectral fitting technique, NBURSTS (Chilingarian et al. 2007a,b), in the basic mode, where a spectrum is approximated by a single simple stellar population (SSP) model. Secondly, we recovered a non-parametric stellar line-of-sight velocity distribution (LOSVD) using an un-broadened stellar population model template from the first step (see method details in Kasparova et al. 2020). The non-parametrical stellar LOSVD reveals a clearly separated double-peaked structure, on which we adjusted a double-Gaussian function to derive kinematics parameters of individual components. We hence spotlight two stellar population components: One is regular, and the second is less prominent but is detected in the region where we see a second gas component as well. We used quality filters based on the uncertainties of the parameters and a sufficient relative flux for the second component to show relevant 2D maps.

The second and fourth columns of Fig. 3 display the stellar velocity and stellar velocity dispersion fields, respectively, for each component. Finally, we re-fitted spectra with NBURSTS, using a two-component approach in which both components have individual and independent kinematics. We considered velocity and dispersion estimates from the previous step as initial guesses. This analysis suggests that both components have similar but slightly different (SSP-equivalent) averaged stellar ages of T₁ ≈ 1.70 Gyr and T₂ ≈ 1.42 Gyr and metallicities of [Fe/H]₁ ≈ −0.24 dex and [Fe/H]₂ ≈ −0.82 dex.

3. Analysis of other data

3.1. Molecular gas content

IRAM 30 m telescope observations were carried out in December 2019 with the Eight MIxer Receiver (EMIR) and the symmetrical wobbler switching mode, in which the secondary mirror switches up to a maximum of ±120 arcsec in azimuth. We operated in dual-sideband mode, and these sidebands were connected to two spectroscopic back ends, WILMA and FTS. J221024.49+114247.0 was integrated for 2.0 h in the E090 band and 2.8 h in the E230 band, enabling the detection of CO(1−0) emission in the former and CO(2−1) emission in the latter, as displayed in Fig. 4. We fitted linear baselines and subtracted them from the spectra. The spectral lines exhibit wide tails and are better adjusted with two Gaussian functions of different widths (see Fig. 4 and the parameters in Table 2) than with a single Gaussian function fit. Firstly, we calculated the CO line luminosity, in K km s⁻¹ pc², as defined in Solomon et al. (1997):

$$\begin{array}{c}\hfill L_{\mathrm{CO}}=3.25\times {10}^7{S}_{\mathrm{CO}}\mathrm{\Delta }V{\nu }_{\mathrm{obs}}^{-2}{D}_{\mathrm{L}}^2{(1+z)}^{-3}\end{array}$$(1)

Fig. 4. Molecular gas. We show the CO(1−0) and CO (2−1) lines observed with the IRAM 30 m telescope in units of main beam temperature, Tmb.

where S_CO = 5 I_CO is the integrated flux in Jy km s⁻¹, I_CO is the line area given in Table 2 in K km s⁻¹ (we considered the sum of both I_CO in order to integrate the whole flux), ν_obs is the line frequency in GHz, and D_L = 422 Mpc is the luminosity distance to the galaxy. Secondly, we derived the mass of molecular gas using M_H2 = α L_CO, where α is the CO-to-H₂ conversion factor, estimated following Freundlich et al. (2019). We hence derived a CO-to-H₂ conversion factor of α = 4 M_⊙/(K km s⁻¹ pc²) and finally found a total molecular gas mass of M_H2 = 3.29 × 10¹⁰ M_⊙, which corresponds to 21% of the total stellar mass. Table 2 gives the properties of the reduced CO spectra.

3.2. Faint structures detected in the r-band MegaPipe/CFHT image

We also used the publicly available r-band stacked image from MegaCam, a wide-field imager mounted on the Canada-France-Hawaii Telescope (CFHT), with very precise astrometric and photometric calibrations, respectively accurate to within 0.15″ and 0.03 mag. By smoothing the image from the MegaCam Image Stacking Pipeline (MegaPipe; Gwyn 2008) using a Gaussian kernel and then subtracting the result from the original image, we spotlight the presence of an obscured dust lane across the galaxy (see the right panel of Fig. 5) and some irregular arms of stars around the core of the system.

Fig. 5. Morphology. Left: composite image (g, r, z bands) from the Legacy survey (Dey et al. 2019). Right: difference between the MegaCam r-band filter image and its smoothed version.

3.3. Radio continuum VLA archives

J221024.49+114247.0 was observed with the Very Large Array (VLA) on 2016 September 30 in A configuration in both frequencies 1.5 GHz and 2.5 GHz (which respectively correspond to the L and the S band), with a 1.3 × 1.3″ resolution for the former and a 0.65 × 0.65″ resolution for the latter (P.I.: J. D. Gelfand). The corresponding fluxes are respectively 8.2 mJy and 3.1 mJy. In parallel, we collected the following integrated fluxes from archival data: 7.4 ± 0.3 mJy at 1.4 GHz from the Faint Images of the Radio Sky at Twenty-Centimeters (FIRST) survey (Becker et al. 1994; Helfand et al. 2015) and 3.1 ± 0.5 mJy at 3 GHz from the Very Large Array Sky Survey (VLASS; Lacy et al. 2020). We derived a synchrotron spectral index of −1.3 and a rest-frame 1.4 GHz luminosity of 1.3 × 10²³ W Hz⁻¹ (Condon et al. 2019). As displayed in the bottom-left panel of Fig. 2, the 1.5 GHz peak coincides with the second Hα component.

4. Results: Properties of the merger

4.1. Kinematics and geometry

4.1.1. Two stellar and gaseous discs

We highlight the superposition of two components along the line of sight in MaNGA data and interpret these spectral features as the two progenitors of an on-going merger event. This interaction involves one galaxy whose gas kinematics appears fairly regular and a second galaxy that is counter-rotating with respect to the first one. As displayed in Fig. 3, the stellar velocity motions for both components are consistent with the velocity fields of the gas. The position angles (PA) of the gas and stellar velocity fields, computed using the algorithms described in Krajnović et al. (2006) (see FIT_KINEMATIC_PA in Cappellari 2017), are similar: PA_gas, 1 = 280.0 ± 0.5° and PA_stars, 1 = 280.0 ± 0.5° for the first component, and PA_gas, 2 = 140.5 ± 2.5° and PA_stars, 2 = 150.5 ± 2.0° for the second component. Using the axes found for the gas kinematics, we estimated lengths of the major and minor axes of, respectively, 41 kpc and 27 kpc for the main galaxy, and of 13 kpc and 9 kpc for the second galaxy.

For the main galaxy, we derived an inclination of i₁ = 49°, assuming a negligible disc thickness (e.g. Hubble 1926) of up to 53° with the correction proposed by Aaronson et al. (1980). For the second component, we derived an inclination i₂ = 46° using the basic formula and i₂ = 50° when applying the correction for disc thickness.

4.1.2. Dust lane

The dust lane displayed in the right panel of Fig. 5 extends from the south-east to the north-west over approximately 11 kpc. We estimated a position angle (PA_dl = 147°) consistent with the position angle derived from the velocity field of the second MaNGA component. As displayed in Fig. 6, this dust lane is compatible with the extinction computed with the Balmer decrement for the second component, with the same inclination, while it is somehow detected in the extinction measured for the first component, supporting the view that the second component lies in front of the system. When compared with the stellar and gaseous component (middle and top panels), it exhibits the same position angle as the velocities but shifted to the north, and hence compatible with an inclination of about 50°.

Fig. 6. Superimposition of the dust lane on the gas velocity maps (top), the stellar velocity maps (middle), and the extinction maps (bottom) derived from Balmer decrement computation.

4.2. Star formation activity

4.2.1. Hα-based SFR

We derive an extinction-corrected star formation rate (SFR) for each component based on the measurements from the double-component decomposition. We corrected the Hα luminosity following Calzetti (2001): L_int(Hα) = L_obs(Hα)10^{0.4k(Hα)E(B − V)}, where L_int is the intrinsic luminosity, L_obs is the observed luminosity, k(Hα) is the reddening curve at the Hα wavelength, and E(B − V) is the dust extinction. To compute this last value, we used the Hα/Hβ ratio and assumed an intrinsic value of 2.86 (corresponding to a temperature of T = 10⁴ K and an electron density n_e = 10² cm⁻³ for a case B recombination; Osterbrock & Ferland 2006). Our derived values based on the MaNGA data, SFR_Hα, 1 = 7.9 M_⊙ yr⁻¹ and SFR_Hα, 2 = 2.3 M_⊙ yr⁻¹, are consistent with the different estimates of the SFR, namely from Brinchmann et al. (2004) and from a spectral energy distribution (SED) fitting based on UV and optical spectra performed by Salim et al. (2016), which agree towards a total value above 10 M_⊙ yr⁻¹.

4.2.2. WISE-based SFR

The SFR derived from Wide-field Infrared Survey Explorer (WISE) data provides 10.57 M_⊙ yr⁻¹ (Salim et al. 2016). It is in excellent agreement with the previous Hα-based estimate, supporting the assumption that there is no significant underestimation of the star formation activity. As discussed later in Sect. 4.2.5, the two galaxies lie on the star formation main sequence.

4.2.3. Radio-based SFR

Following Condon et al. (2019) and assuming the radio flux is due to star formation only, we also derived a radio-based SFR that could be as large as 160 M_⊙ yr⁻¹. As mentioned in Sect. 3.3, the slope computed from the four flux densities (FIRST, both bands of VLA, and VLASS) is very steep. Although this value of −1.3 is not sufficient to assess whether the emission originates from AGNs or star formation, it is typical of a synchrotron emission rather than a free-free one. The discrepancy between the optical-based and the radio-based SFR estimates is due to either a severe underestimation of the extinction in the optical estimate or a hidden AGN or Low-Ionisation Nuclear Emission-line Region (LINER) contributing to the radio flux. However, the good agreement between the optical Hα-based SFR and the infrared WISE-based SFR discussed above rules out the first hypothesis.

As displayed in the bottom-left panel of Fig. 2, the peak at 1.5 GHz (as well as the peak at 2.5 GHz) is coincident with the second component Hα flux. Given the large SFR observed for the main galaxy, the radio flux due to star formation would be expected to peak on this galaxy. The fact that the radio flux is centred on the second component is a strong argument in favour of a hidden AGN and supports the notion that most of this radio flux (more than 90%) is not due to star formation.

4.2.4. Off-centred star formation

For the main galaxy, the peak of star formation detected in Hα is off-centred with respect to the kinematic centre. The velocity dispersion peaks at the centre of the rotation pattern, but some large gas velocity dispersions are detected in the north-west part of the disc. Indeed, while the stellar velocities are in good agreement with gas velocities, the stellar velocity dispersions are shifted with respect to those observed for the gas. Moreover, we cannot exclude that the extinction correction of Hα emission is not sufficient to trace the actual ionised gas emission in the centre of the main galaxy. We might thus miss star formation at its centre. The stellar velocity dispersion might also be affected by an important amount of dust in the centre, which would explain why we detect high values in the north-west region of the galaxy where there is less dust. The presence of clumps and filament-like structures, observed in the dust extinction maps (Fig. 2) and in the velocity dispersion map of the primary galaxy (Fig. 3), also suggests the presence of underlying velocity gradients due to recent gas accretion triggered by the interaction.

4.2.5. Two main-sequence galaxies

Figure 7 summarises the main star-forming properties of the two galaxies. They have a regular star-forming activity as both lie on the main sequence. For comparison, we display the contours based on galaxies taken from the RCSED (Chilingarian et al. 2017). This catalogue encompasses 800 299 low- and intermediate-redshift galaxies (0.007 < z < 0.6) and provides multi-wavelength spectrophotometric measurements from various surveys (GALEX, SDSS, and UKIDSS) as well as value-added data (k corrections, star and gas kinematics, and extinction estimates) for each object. Globally, the two detected components currently do not exhibit any excess of star formation. This result is compatible with the molecular gas detection discussed in Sect. 4.5. Lastly, the excitation diagram explored next also supports the star-forming and composite excitations for most regions of these two galaxies.

Fig. 7. SFR vs. stellar mass plot. The black star corresponds to the total system, relying on the values provided by Salim et al. (2016) for J221024.49+114247.0. The blue and red stars refer to the main and secondary component, respectively, based on our SFR and stellar mass computations. The contours correspond to representative galaxies of the RCSED sample (Chilingarian et al. 2017), classified using Domínguez Sánchez et al. (2018).

4.3. BPT excitation diagram

When considering the Baldwin, Phillips & Telervich diagnostic diagram (introduced in Baldwin et al. 1981) for each peak separately (cf. Fig. 8), most of the spaxels lie in the composite and star-forming regions for both peaks, but the second peak displays a systematically higher [OIII]/Hβ ratio than the first one. This supports our two-component analysis and is in line with the fact that the redshifted component has a smaller mass (and lower metallicity; e.g. Kewley et al. 2019; Curti et al. 2020). While the two discs are mostly star forming, the BPT maps (Fig. 8) display very localised LINER-type emission.

Fig. 8. Ionisation properties of the gas for the two components of J221024.49+114247.0. Left and middle panels: BPT (Baldwin et al. 1981) maps for both fitted components. Right panel: BPT diagram showing the position of the spaxels. The S/N threshold is set at 3. The spaxels where the S/N of at least one emission line is below the corresponding threshold are masked. Moreover, for the redshifted component, the masks described in Sect. 2.1 are applied. The last panel displays the LINER, AGN, composite, and star-forming excitation regions as defined by Kewley et al. (2001), Kauffmann et al. (2003), and Schawinski et al. (2007). The grey contours refer to galaxies representative of the RCSED sample.

Based on the [OIII]λ5008 emission, it appears that the double-peaked profile for this line is reliably detected in a wider part of the MaNGA field of view. The spaxels identified as the secondary component do not seem to be associated with the main velocity gradient (their velocity is below 0 km s⁻¹) and show a higher excitation compared to the other spaxels related to this galaxy. This is reflected in their position on the BPT diagram, slightly shifted with respect to the main distribution of the other points and clearly belonging to the composite region (or even to the LINER region for three of them). The integrated [OIII]λ5008 luminosities for each component within an area equivalent to the PSF (a circular region of about 2.5″) around the respective emission peak are L_[OIII],1 = 1.1 × 10⁷ L_⊙ and L_[OIII],2 = 8.3 × 10⁶ L_⊙. Even though we applied an extinction correction to derive these values, this might not be sufficient to trace the total [OIII]λ5008 emission in these inner parts of the galaxies. The values we find are higher than 10^6.5 L_⊙, the limit between weak and powerful AGNs considered by Alonso et al. (2007), which is consistent with their conclusion that most of the AGNs in pairs show an [OIII] luminosity above this threshold.

In Fig. 9 we also study stacked histograms of [OIII] and Hβ fluxes. We both stacked the spectra of all spaxels and stacked the models of the blue and red components. Interestingly, we find some high velocity residuals that escaped the adjustment. They typically correspond to expected velocity fields in the outskirts of both discs, where the S/N is weak.

Fig. 9. Stacked ionised-gas lines. Top and bottom panels: stacked Hβ and [OIII]λ5008 emission lines, respectively, from all MaNGA spaxels in black. We show with blue and red the stacked blueshifted and redshifted components of our fitted model, only including spaxels with detected lines. Below the stacked emission lines, we show the residuals that exhibit excesses at large velocities.

4.4. Possible evidence of a radio jet

Above we stated that the radio flux seems to be associated with the smaller galaxy. The radio emission from VLA observations is extended: We computed an integrated flux, called F_int hereafter, by summing up the emission over an area of 30 × 30 pixels around the flux peak and divided the derived value by the beam surface. We finally find F_int/F_peak = 8.7 (respectively 5.1) at 2.5 GHz (respectively 1.5 GHz) with a 0.65″ (respectively 1.3″) resolution. Given the radio flux excess with respect to the expected SFR, we could argue that this is a weak radio jet.

4.5. Molecular gas content

The single-dish CO observations reveal an important reservoir of molecular gas within the system, both in CO(1−0) and CO(2−1). The inferred total H₂ mass represents 21% of the stellar mass of the system. The Kennicutt–Schmidt relation shows that the system is located on the main sequence.

When we discriminate the flux below 0 km s⁻¹ from the one above this value, as for the optical data, we find that the galaxy detected in the negative velocities has a higher CO(1−0)/CO(2−1) ratio than the secondary component. The beam sizes are different for both transitions, but since the CO(2−1) beam covers ∼23 kpc, we consider that we do not miss any significant emission outside the beam coverage. Thus, we can associate this ratio discrepancy with different excitations.

4.6. Gas metallicity and post-starburst signature

The maps of the gas metallicity are shown in the right panels of Fig. 2. The smaller component appears to have a lower gas-phase oxygen abundance than the main galaxy, but this is not so clear when one computes averaged metallicities 12 + log(O/H)₁ = 8.65 and 12 + log(O/H)₂ = 8.58. However, the spatial distribution is very different. While the small component has a metallicity gradient roughly centered on the kinematic centre, the main galaxy exhibits a large gas metallicity region east of the star formation peak. This disordered metallicity gradient can be related to the gas accretion from the low-metallicity. The main galaxy follows the mass-metallicity relation when we consider this off-centred area, while it exhibits a lower-than-expected metallicity in the central parts, as explained in Sect. 4.7.4 and shown in Fig. 10. Moreover, the galaxy spectra in these off-centred regions exhibit a significant Balmer absorption line that is attributed to a stellar population of ∼1 Gyr (e.g. Pawlik et al. 2018), typical of post-starburst regions. On the contrary, the spectra extracted from the regions where the computed gas-phase oxygen abundance is below the solar metallicity do not show such an important Balmer absorption line. Typical spectra extracted from these two different types of regions are displayed in Fig. B.1. Unlike in Goto (2007), the galaxies are still forming stars. Also, the starburst is quite old (1 Gyr) compared to the systems observed in, for example, Poggianti & Wu (2000).

Fig. 10. Mass-metallicity relations for different values of total SFR, as defined in Curti et al. (2020). The SFRs of both components are the ones computed from our two-component analysis and given in Sect. 4.2.1. The markers are colour-coded with the SFR of the corresponding host galaxies. For the primary galaxy, different 3 arcsec circular apertures have been considered, as discussed in Sect. 4.7.4.

However, we do not see explicit or strong signs of a recently formed stellar population in the optical spectrum, but the S/N of the data makes it difficult to provide definitive conclusions based on individual bins. Ultraviolet data with a good spatial resolution would help us conduct a qualitative analysis of these regions.

The irregular metallicity distribution can be compared with other merging systems, such as NGC 4038/4039. This system was studied by Gunawardhana et al. (2020), who detected star-forming regions with gas more enriched than the rest of the galaxy.

4.7. Mass ratio estimates

4.7.1. Stellar mass ratio

To determine the light fraction between both detected stellar components, we performed a modelling of the stellar LOSVD. The recovered LOSVD (see Sect. 2.2) is approximated in each spatial bin by two Gaussian components. Each of these components is represented by a thin inclined exponential disc. We used arbitrary initial guesses for the different characteristics of these profiles and eventually found an optimal set of parameters to reproduce the modelled LOSVD. We obtain a light fraction of 6.1 between the main component and the secondary one. Using stellar ages and metallicities derived with NBURSTS, as mentioned in Sect. 2.2, we determined the mass-to-light (M/L) ratios for the r-band filter: M/L_r, 1 = 0.63 and M/L_r, 2 = 0.43. Given that the light fraction is 6.1, the stellar mass fraction is as follows: 6.1 (M/L_r, 1)/(M/L_r, 2) = 9. This mass ratio enabled us to infer stellar masses for each of the detected component: According to Salim et al. (2016), the total stellar mass of the system is 1.59 × 10¹¹ M_⊙, which, when considering the stellar mass ratio described above, gives: M_⋆, 1 = 1.43 × 10¹¹ M_⊙ and M_⋆, 2 = 1.59 × 10¹⁰ M_⊙.

4.7.2. Dynamical mass ratio

For each galaxy, the inclination angle computed in Sect. 4.1 can be used to calculate the respective rotation velocity as in Cortese et al. (2014) and Aquino-Ortíz et al. (2018) following V_rot = W/2 sin(i), where W is defined as the difference between the 90th and 10th percentile points of the velocity histogram. We derive a rotation velocity of V_rot, 1 = 236 km s⁻¹ for the main galaxy and V_rot, 2 = 140 km s⁻¹ for the secondary one. Then we computed the dynamical mass corresponding to the first MaNGA component as follows: $M_{\text{dyn},\text{1}}=V_{\text{rot},\text{1}}^2{R}_1/G$, where R₁ is the maximal radius where Hα emission is detected. We find M_dyn, 1 = 2.65 × 10¹¹ M_⊙. Similarly, the dynamical mass of the second component is defined as: $M_{\text{dyn},\text{2}}=V_{\text{rot},\text{2}}^2{R}_2/G$. We derive M_dyn, 2 = 2.93 × 10¹⁰ M_⊙. These results give a mass ratio of M_dyn, 1/M_dyn, 2 ∼ 9.

4.7.3. Uncertainties

To evaluate the robustness of the previous estimate of the merger mass ratio, we varied the parameters used in the computations of the stellar and the dynamical masses. As the interaction between both objects certainly perturbs both galaxies, we can consider that it is relevant to include the contribution of the velocity dispersion σ_gas in the dynamical mass estimate. Based on the S_0.5 parameter defined and discussed in Aquino-Ortíz et al. (2018), we calculated the dynamical mass using the proxy they introduce as $M_{\text{dyn}}=\eta R{S}_{0.5}^2/G$. The inferred dynamical masses are M_dyn, 1 = 6.98 × 10¹¹ M_⊙ and M_dyn, 2 = 5.80 × 10¹⁰ M_⊙, corresponding to a mass ratio M_dyn, 1/M_dyn, 2 of 12. In parallel, for the stellar masses, we can force both inclinations of the modelled exponential disc profiles based on the inclination angles computed using the gas velocity fields. This leads to a light fraction of 7.3. Relying on the mass-to-light ratios given in Sect. 4.7.1, the corresponding stellar mass ratio is ∼10. Lastly, one can note that considering a large inclination angle for the secondary galaxy provides a dynamical mass similar to the stellar mass, constraining an inclination smaller than 80°.

Finally, all the approaches converge on a mass ratio within the visible radius of about 9:1 for MaNGA 1−114955.

4.7.4. Mass-metallicity relation

In order to determine if the metallicities of these two galaxies, which both lie on the star formation main sequence, are perturbed, we tried to position them on the stellar mass-SFR-gas metallicity relation proposed by Curti et al. (2020) (based on Mannucci et al. 2010). Figure 10 shows the mass-metallicity relations for different SFR bins as described in Curti et al. (2020). It should be noted that the gas metallicity computed here corresponds to the gas metallicity computed in the SDSS 3 arcsec (central) fibre. We used the calibration defined in this work to compute the metallicity, based on the O₃N₂ calibrator, within circular regions with a diameter of 3 arcsec. On the one hand, the second galaxy lies almost on the curve parametrised for objects with a similar SFR. On the other hand, the properties of the first galaxy are consistent with the relation only if we compute the metallicity over the off-centred region with high gas-phase oxygen abundance values (see the top-right panel of Fig. 2). If we compute the metallicity around the kinematic centre or the star formation peak, the value is too low compared to the mass of the galaxy and thus does not follow the relation. We can argue that this primary galaxy has probably accreted metal-poor gas in its centre through the interaction.

5. Discussion

5.1. Pre-coalescence merger

We identify J221024.49+114247.0 as a merging system caught in action. We observe two distinct counter-rotating discs, aligned along the line of sight, with a line-of-sight velocity difference of ΔV = 450 km s⁻¹. We find one main galaxy with an inclination of i₁ = 49° and a second galaxy of i₂ = 46° as a smaller counterpart. Using high-resolution imaging from the Megacam-CFHT footprint, we identified a dust lane that is aligned with the stellar component of the smaller galaxy and is associated with measured extinction in the main galaxy, suggesting a scenario where the small galaxy is positioned in front of the main galaxy along the line of sight, as shown in Fig. 11. These two disc components are both detected in the ionised gas and in the stellar continuum emission, and they have not yet been destroyed by the interaction. We estimated a mass ratio of 9:1 using the stellar masses and dynamical masses. For both galaxies, we find off-centred peaks in stellar velocity dispersion, reminiscent of possible accretion onto the disc and an underlying filamentary structure.

Fig. 11. Schematic illustration of the revealed merging system. Top: configuration where the near side of the secondary galaxy is on the north side. Bottom: configuration where the near side is on the south side of the major axis. Left panels: discs as projected on the sky. Right panels: a ‘3D view’, with a different viewing angle, for which discs can be represented with no superposition. Top-left panel: highlighting of a potential dust lane in Fig. 5.

We show schematic representations of two possible configurations of the observed objects in Fig. 11: Each row of the figure shows one of these configurations, viewed as projected on the sky (left panels) or with a viewing angle (the same for the two rows) for which the discs can be represented with no superposition (right panels, ‘3D views’). The line of sight along which the system is observed with MaNGA is parallel to the (Oy) axis (the MaNGA observer is symbolised with the eye to the right). The disc on the left of the 3D views represents the primary galaxy, and the one on the right, the secondary. The rotation of each disc is shown with a spin (black arrow) and corresponding schematic trailing spiral arms. The coloured arrows represent the line-of-sight velocity for the MaNGA observer (blue- or red-shifted with respect to each kinematic centre), and the MaNGA velocity maps are schematically represented on the (Oxz) plane of the right panels (with different colours corresponding to each kinematic centre, as in Fig. 3). We note that the distance between the two galaxies is unknown (hence the break in the line of sight). For a given disc inclination, part of one side of the major axis of the projection of the disc on the sky is either farther away than this major axis or closer to it (near side). We can determine the near side of the primary galaxy with the help of its spiral arms, seen in Fig. 5. We assumed that these spiral arms are trailing, which gives the sense of rotation on the sky. Combined with the line-of-sight velocities, this determines the inclination of the disc plane, with a near side north-east of the major axis. The near side of the secondary galaxy is harder to determine. The presence of the dust lane parallel to the kinematic axis but shifted to the north could support an argument that the near side is on the north side, as in the top panel of Fig. 11. Alternatively, the bottom panel shows a configuration compatible with the observation of an ionisation cone north-east of the major axis of the secondary disc, which would mean that the near side is on the south side of the major axis.

5.2. Possible outflows driven by central activity

The interaction might also have triggered gas accretion on the supermassive black hole (SMBH) and a hidden AGN detected in radio. This is compatible with the discussion in Shah et al. (2020) regarding the possibility of low-level AGN enhancement in galaxy interactions. Using BPT diagrams, we find that the two galaxies are located in different regions, associated with different states of ionisation. On the one hand, the star-forming and composite excitation occupies different regions, probably due to the different galaxy masses, while the second galaxy shows some regions with LINER excitation, suggesting a different mechanism that is potentially linked to the radio flux found to be associated with this galaxy. This might indicate that the merger process can induce different mechanisms, shifting the two components on the BPT diagram, as discussed in Maschmann & Melchior (2019).

The molecular gas spectra obtained for both galaxies exhibit a large range of velocities of about 1000 km s⁻¹. By comparing these spectra to the stacked Hβ and [OIII] emission lines from the optical spectra in Fig. 12, we observe that the wide wings of CO emission coincide with the wings observed in the ionised gas. This accounts for the rotation patterns of the two discs. In parallel, we detect, in the velocity range −100 to 0 km s⁻¹, a relative excess of molecular gas corresponding to the peak of star formation activities. To understand how the molecular gas behaves in the vicinity of the possible AGN identified in radio, higher-resolution observations of molecular gas with the NOrthern Extended Millimeter Array (NOEMA) are needed.

Fig. 12. Re-scaled CO spectra and fits of the stacked Hβ and [OIII] lines that display a similar range of velocities.

The wide range of velocities detected in both molecular and ionised gas might correspond to a gaseous outflow, either due to the star formation triggered by the interaction or to the feedback due to the AGN in the small galaxy. It has been shown that radio jets, even when compact, can create major turmoil in galaxies, enhancing the velocity width (Venturi et al. 2021). Figure 9 reveals that the wide wings are not detected in individual spectra, but appear only by stacking. This turbulent or outflowing gas must be spread over some diffuse regions and only at low intensity.

5.3. Evidence of two pericentre passages

The total SFR of the entire system is estimated to be 13.3 M_⊙ yr⁻¹ (Salim et al. 2016). Radio observations suggest a much higher SFR of more than 80 M_⊙ yr⁻¹, which is probably biased by a contribution of an obscured AGN. Using the extinction-corrected Hα emission line, we located a concentration of ongoing star formation in the main galaxy, which is off-centred with respect to the kinematic centre by about 5 kpc in projection. The smaller galaxy exhibits a strong peak in star formation, located at its kinematic centre, and is associated with an extended radio emission.

It seems that the external gas of the smaller companion has been stripped by the main galaxy, which has probably accreted this gas onto its disc. It is hence probable that the progenitor of the smaller galaxy was originally more massive. The ignition of the AGN, somehow related to the gas present in the central starburst, is due either to accretion of gas (from the primary) or to the original central distribution of gas, and is not necessary related to the harassment. We also observe a peak in the gas-phase metallicity of the main galaxy located to the east of its main star formation region. This high-metallicity region, located more than 6 kpc away from the kinematic centre, is associated with significant Balmer absorption lines imprinted in the galaxy spectrum (see Fig. B.1). This feature could trace past starburst events not older than 1 Gyr (e.g. Pawlik et al. 2018). Indeed, this clearly asymmetric metallicity distribution is not expected in a relaxed system. We argue that this could be the signature of a first pericentre passage, while the current star formation (in the north-west part of the main galaxy) could be associated with a recent pericentre passage. Star formation induced by the interaction of galaxies varies significantly (in terms of highest SFR, duration of high SFR periods, and integrated star formation on the duration of the merger) as a function of the orbital parameters and the inclination of the discs of the progenitors (e.g. Di Matteo et al. 2008). It is thus difficult to infer the geometry of the interaction (prograde or retrograde orientation of the spins of the galaxies) from the observed star formation.

5.4. Comparison with other merging systems

The interacting galaxies presented in this work can be compared to the less massive system Mrk 739, which hosts a double AGN (Tubín et al. 2021). In this case, two AGNs have been identified in X-ray (Koss et al. 2011). Interestingly, in this source the ionised gas distribution clearly peaks at the position of the two AGNs. The gas is concentrated in the centre, which explains the ignition of both AGNs. As the system is more face-on with a projected nuclear separation of 3.4 kpc, it is not obvious if there are two rotating distinct discs, such as those we observe for J221024.49+114247.0. The authors noted that the progenitors were a star-forming galaxy and an old elliptical.

In our case, the main galaxy has a larger SFR than its companion, but this is probably triggered by accretion of the stripped gas. Interestingly, this gas did not reach the galaxy centre and was not able to fuel the SMBH. Using a simulation of NGC 2623, Prieto et al. (2021) explained that the first passage can trigger extended star formation, such as we observed here. He et al. (2020) discussed the two galaxies in Arp 240 with a projected distance of about 40 kpc. They claim that the low depletion time (100 Myr) is an argument for a late-stage merger. In the case of J221024.49+114247.0, the global depletion time (3.3 Gyr) is typical of normal star-forming galaxies (e.g. Saintonge et al. 2011).

6. Conclusions

We have presented a multi-wavelength analysis of a pre-coalescence merger, J221024.49+114247.0, at z = 0.09. Originally identified with a double-peaked feature in the SDSS fibre and a disrupted merger morphology, this galaxy has been observed with MaNGA. We developed a fitting procedure based on a double-Gaussian function and applied it throughout the MaNGA field of view. We have been able to highlight the presence of two counter-rotating discs along the line of sight in gas as well as in stars. The difference of their systemic velocity is 450 km s⁻¹, and they are aligned along the line of sight. A modelling of the stellar LOSVD enables us to derive a stellar mass ratio between the main and the secondary component of ∼9, relying on the mass-to-light ratio estimates. The mass ratio based on dynamical mass estimates is in good agreement with that derived from the stellar distribution analysis.

The main galaxy exhibits off-centred regions of star formation. A large region (at ∼6 kpc from the kinematic centre) exhibits a large gas metallicity and a Balmer absorption line in the galaxy spectra, possibly reminiscent of a 1 Gyr old starburst, while the current SFR detected in Hα is also off-centred with respect to the kinematic centre. We argue that this is due to two different pericentric passages. This more massive galaxy is probably more gas rich than its companion, but the gas has not collapsed in the centre. The minor galaxy exhibits an excess of extended radio flux centred on its kinematic centre with respect to the expected SFR. Its kinematic centre corresponds to a standard velocity dispersion and metallicity gradient. We argue that its external gas has been stripped. Large molecular gas velocities are compatible with the ionised gas rotation velocities, and we cannot conclude on the possibility of an outflow.

We propose to extend this work with NOEMA/IRAM interferometric observations in order to study the molecular gas distribution and the sub-millimetre continuum distribution.

Acknowledgments

We are most grateful to the IRAM-30 m team who supported us for these observations. IK acknowledges the support from the Russian Scientific Foundation grant 19-12-00281 and the Interdisciplinary Scientific and Educational School of Moscow University “Fundamental and Applied Space Research”. We thank the anonymous referee for constructive comments. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the US 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 website 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, the Korean Participation Group, 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. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.

References

Appendix A: Gas spatial properties obtained through our fitting procedure

A.1. Flux maps

Fig. A.1. Respective Hβ, [OIII]λ5007, Hα, and [NII]λ6585 flux maps for both fitted components. The S/N threshold is set at 3 (spaxels with S/N below this value are masked). Moreover, for the second component, all the spaxels that do not meet the criteria defined in Sect. 2.1 are masked. The maps shown in the right column are the respective fluxes but corrected for extinction.

Appendix B: Spectra from regions with different extreme metallicities

Fig. B.1. Different extracted spectra and their corresponding fits. Left panel: spectrum of a bin with a high gas-phase oxygen abundance, and right panel: a bin where the computed gas-phase oxygen abundance is low. Left panel: an important Balmer absorption line below Hβ can be seen.

All Tables

All Figures

	Fig. 1. Velocity dispersion of the Hβ emission line derived by the MaNGA DAP. The MaNGA field of view is represented as a grey footprint. The included panel shows a portion of the continuum-subtracted spectrum from one spaxel where the model (in orange) does not manage to fit the spectral features.
In the text

Fig. 2. Gas properties derived from our multi-component approach. First and second rows: maps for the first and second components detected in MaNGA data, respectively. First column: Hα extinction-corrected flux (in erg s−1 Å−1 cm−2 per spaxel), second column: extinction computed from the Balmer decrement, and last column: oxygen gas-phase abundance derived using the O3N2 calibrator. The black crosses indicate the position of the extinction-corrected Hα flux peak for the represented component. The MaNGA PSF is displayed as a hatched grey circle in the bottom-left corner of the panels.

In the text

Fig. 3. Gas and star kinematics. First row: maps for the first component, and second row: second component detected in MaNGA data. First and second columns: respective ionised gas and stellar velocity fields (in km s−1). The dotted black lines refer to the computed position angles. Third and fourth columns: respective velocity dispersion fields (in km s−1) for the gas and the stars. As in Fig. 2, the black crosses indicate the position of the extinction-corrected Hα flux peak for the represented component.

In the text

	Fig. 4. Molecular gas. We show the CO(1−0) and CO (2−1) lines observed with the IRAM 30 m telescope in units of main beam temperature, Tmb.
In the text

	Fig. 5. Morphology. Left: composite image (g, r, z bands) from the Legacy survey (Dey et al. 2019). Right: difference between the MegaCam r-band filter image and its smoothed version.
In the text

	Fig. 6. Superimposition of the dust lane on the gas velocity maps (top), the stellar velocity maps (middle), and the extinction maps (bottom) derived from Balmer decrement computation.
In the text

Fig. 7. SFR vs. stellar mass plot. The black star corresponds to the total system, relying on the values provided by Salim et al. (2016) for J221024.49+114247.0. The blue and red stars refer to the main and secondary component, respectively, based on our SFR and stellar mass computations. The contours correspond to representative galaxies of the RCSED sample (Chilingarian et al. 2017), classified using Domínguez Sánchez et al. (2018).

In the text

Fig. 8. Ionisation properties of the gas for the two components of J221024.49+114247.0. Left and middle panels: BPT (Baldwin et al. 1981) maps for both fitted components. Right panel: BPT diagram showing the position of the spaxels. The S/N threshold is set at 3. The spaxels where the S/N of at least one emission line is below the corresponding threshold are masked. Moreover, for the redshifted component, the masks described in Sect. 2.1 are applied. The last panel displays the LINER, AGN, composite, and star-forming excitation regions as defined by Kewley et al. (2001), Kauffmann et al. (2003), and Schawinski et al. (2007). The grey contours refer to galaxies representative of the RCSED sample.

In the text

	Fig. 9. Stacked ionised-gas lines. Top and bottom panels: stacked Hβ and [OIII]λ5008 emission lines, respectively, from all MaNGA spaxels in black. We show with blue and red the stacked blueshifted and redshifted components of our fitted model, only including spaxels with detected lines. Below the stacked emission lines, we show the residuals that exhibit excesses at large velocities.
In the text

Fig. 10. Mass-metallicity relations for different values of total SFR, as defined in Curti et al. (2020). The SFRs of both components are the ones computed from our two-component analysis and given in Sect. 4.2.1. The markers are colour-coded with the SFR of the corresponding host galaxies. For the primary galaxy, different 3 arcsec circular apertures have been considered, as discussed in Sect. 4.7.4.

In the text

Fig. 11. Schematic illustration of the revealed merging system. Top: configuration where the near side of the secondary galaxy is on the north side. Bottom: configuration where the near side is on the south side of the major axis. Left panels: discs as projected on the sky. Right panels: a ‘3D view’, with a different viewing angle, for which discs can be represented with no superposition. Top-left panel: highlighting of a potential dust lane in Fig. 5.

In the text

	Fig. 12. Re-scaled CO spectra and fits of the stacked Hβ and [OIII] lines that display a similar range of velocities.
In the text

	Fig. A.1. Respective Hβ, [OIII]λ5007, Hα, and [NII]λ6585 flux maps for both fitted components. The S/N threshold is set at 3 (spaxels with S/N below this value are masked). Moreover, for the second component, all the spaxels that do not meet the criteria defined in Sect. 2.1 are masked. The maps shown in the right column are the respective fluxes but corrected for extinction.
In the text

	Fig. B.1. Different extracted spectra and their corresponding fits. Left panel: spectrum of a bin with a high gas-phase oxygen abundance, and right panel: a bin where the computed gas-phase oxygen abundance is low. Left panel: an important Balmer absorption line below Hβ can be seen.
In the text