© ESO 2019

1. Introduction

Powerful active galactic nuclei (AGN) and starburst (SB) -driven outflows are routinely invoked to link the growth of supermassive black holes (SMBH) and galaxies: these phenomena, which deposit energy and momentum into the surrounding medium, are expected to affect the physical and dynamical conditions of infalling gas and thus regulate the formation of new stars and the accretion onto the SMBH (see Somerville & Davé 2015 for a recent review). In particular, feedback phenomena are thought play an important role in the so-called starburst-quasar (SB-QSO) evolutionary sequence framework (e.g. Sanders et al. 1988), where a dusty, SB-dominated system, arising from a merger, evolves into an optically luminous QSO (e.g. Hopkins et al. 2008; Menci et al. 2008).

Most star-forming galaxies and AGN show evidence of powerful galaxy-wide outflows. Spatially resolved optical, infrared (IR), and millimenter (mm) spectroscopic studies are revealing the extension, morphology, and energetics related to the ejected multi-phase material, in both the local Universe (e.g. Arribas et al. 2014; Contursi et al. 2013; Feruglio et al. 2015; Harrison et al. 2014; Revalski et al. 2018; Sturm et al. 2011; Venturi et al. 2018a) and at high redshift (e.g. Brusa et al. 2016, 2018; Carniani et al. 2016, 2017; Cresci et al. 2015a; Föerster-Schreiber et al. 2014; Genzel et al. 2014; Harrison et al. 2012, 2016; Newman et al. 2012; Perna et al. 2015a,b, 2018a; Vietri et al. 2018). The study of the different phases of the outflowing gas in individual galaxies, which is essential to constrain feedback model prescriptions (see e.g. Cicone et al. 2018; Harrison et al. 2018), is limited to a very small number of targets, however (e.g. Brusa et al. 2015a, 2018; Carniani et al. 2017; Feruglio et al. 2015, 2017; Tombesi et al. 2015; Veilleux et al. 2017).

While most of the evidence of AGN outflows in nearby galaxies comes from the detection of blue-shifted (and red-shifted) outflow components of ionised and molecular gas (generally traced by [O III]λ5007 and CO emission lines), several studies have indicated that only ∼1% of nearby AGN show evidence of neutral outflows, which is usually traced by resonant sodium lines at 5890 and 5896 Å (Perna et al. 2017a and references therein). Neutral outflows are instead easily detected in ultra-luminous infrared galaxies (ULIRGs) and sources exhibiting concomitant star formation (SF) and AGN activities (Nedelchev et al. 2017; Cazzoli et al. 2016; Rupke et al. 2005a, 2017; Sarzi et al. 2016). These findings could suggest that AGN activity has no impact on the neutral gas (e.g. Bae & Woo 2018; Concas et al. 2017), or alternatively, that neutral outflows can be observed only in obscured AGN in which large reservoirs of cold gas are still present and the outflow processes can affect all the gas phases (Perna et al. 2017a).

In order to distinguish between the two different scenarios, high-quality spatially resolved spectroscopic observations are required. These early results, which are generally based on spatially integrated data, could have been affected by technical and observational limitations: (i) spectra with high signal-to-noise ratios (S/N) are required to trace Na ID kinematics, (ii) a detailed analysis is needed to distinguish between contributions from the interstellar medium (ISM) and stars to the sodium lines, (iii) AGN continuum emission can easily over-shine the Na ID features, and (iv) weak neutral outflows can remain undetected in spatially integrated spectra. All these issues could suggest that the incidence of AGN-driven neutral outflows has been so far underestimated, as also indicated in Rupke & Veilleux (2013, 2015) and Rupke et al. (2017). The ubiquitous presence of molecular outflows in AGN systems (e.g. Fiore et al. 2017 for a recent compilation) indeed indicates that AGN winds can successfully accelerate the cold (T = 10² − 10³ K) ISM components.

In Perna et al. (2017a,b) we presented the spectroscopic analysis of 3″ integrated Sloan Digital Sky Survey (SDSS) spectra of ∼650 X-ray detected AGN. We derived the incidence of ionised outflows (∼40%) and studied the relations between outflow velocity and different AGN power tracers (SMBH mass, [O III]λ5007, and X-ray luminosity). We also derived the incidence of neutral outflows by studying the sodium Na IDλλ5890,5896 transitions. In particular, the north-west galaxy in the Mkn 848 system, a binary merger at z ∼ 0.04, is the only source for which the concomitant presence of outflowing gas in the ionised ([O III]λ5007) and neutral phases has been revealed. Our finding is therefore consistent with previous works, which showed a neutral outflow incidence ≲1% in spatially integrated spectra of local AGN.

Mkn 848 represents an ideal laboratory for feedback studies, as simulations show that gas-rich galaxy mergers can provide the conditions for intense SF and AGN activity, which is essential for the production of powerful outflows. The two merging galaxies are above the so-called main sequence (MS) of normal star-forming galaxies (e.g. Whitaker et al. 2012), with star formation rates (SFR) of ∼80 M_⊙ yr⁻¹ and ∼10 M_⊙ yr⁻¹ for the north-west and south-east galaxies, respectively (Yuan et al. 2018), to be compared with an SFR_MS ∼ 2 M_⊙ yr⁻¹ of local galaxies with similar stellar masses. Moreover, their optical spectra are associated with starburst-AGN composite emission lines (e.g. Yuan et al. 2010, using standard Baldwin–Phillips–Terlevich (BPT) diagnostics; Baldwin et al. 1981). However, such line ratios, when associated with ULIRGs or interacting galaxies, might also be due to extended shocks, as suggested by Rich et al. (2015). Even if the presence of active SMBHs cannot be definitely confirmed by simple optical line ratio diagnostics, Mkn 848 still represents an optimal target for studying the SB-QSO paradigm, as it allows us to test the extreme conditions of the ISM during a merging event.

In this paper, we present a kinematic and physical analysis of Mkn 848 obtained from MaNGA (Mapping Nearby Galaxies at APO) observations, which have been publicly released as part of the SDSS-IV Data Release (Blanton et al. 2017). The manuscript is organised as follows. In Sect. 2 we present the ancillary target data and discuss the properties of Mkn 848 that we derive from available multi-wavelength information. In Sect. 3 we show the MaNGA data analysis, while Sect. 4 displays the spatially resolved spectroscopic results. In Sects. 5 and 6 we estimate extinction and hydrogen column density. In Sect. 7 we derive the plasma properties related to systemic and outflowing gas, while in Sect. 8 we investigate the dominant ionisation source for the emitting gas. Finally, we present the outflow energetics in Sect. 9 and summarise our results in the Conclusion. We adopt the cosmological parameters H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7. The spatial scale of 1″ corresponds to a physical scale of ∼0.8 kpc at the redshift of Mkn 848.

2. Ancillary data for Mkn 848

The source Mkn 848 consists of two merging galaxies with similar stellar masses (log(M_*/M_⊙)≈10.4; Yuan et al. 2018). In Fig. 1 we show the optical, near-infrared, and three-colour images of the system from the Hubble Space Telescope (HST), as well as a Chandra X-ray cutout. The two nuclei have a projected separation of 7.5 kpc. Two long, highly curved tidal tails of gas and stars emerge from the north-west (NW) and south-east (SE) galaxies. One arm (T1 in Fig. 1), curving clockwise, stretches from the top of the image and interlocks with the other arm (T2), which curves up counter-clockwise from below. The gravitational interaction between the galaxies may also be responsible for the shell-like region in between the two nuclei (S in the figure). The three-colour image also shows an extended dust lane in the SE galaxy and strong absorption in the NW and SE nuclear regions.

Fig. 1. Chandra 0.5–8 keV cutout (left panel) and HST images with colours in log scale. In the Chandra image, blue contours represent F435W fluxes. In the F435W image, we labelled the two major tidal tails, T1 and T2, that extend over tens of kpc and possibly overlap on the right of the SE nucleus in the figure. A shell-like region (S) is also indicated. In the three-colour cutout, the white boxes show the GMOS FOV from which previous spatially resolved kinematic results have been obtained by Rupke & Veilleux (2013). The white contour represents the S/N = 40 level in the stellar continuum map at ∼5550 Å and is taken as a reference in the next figures. The MaNGA field of view is also shown with white dashed lines.

Available multi-wavelength information from the UV to the IR has recently been analysed by Yuan et al. (2018): they used photometric and spectroscopic (MaNGA) data to examine the dust attenuation and the star formation activity in six different regions that cover the nuclear cores and tidal arms. The analysis revealed a very recent (∼100 Myr) interaction-induced starburst in the central regions (especially in the NW nucleus), but also signs of earlier (∼500 Myr) star formation activity in the tidal area. The authors also found that the AGN fraction to the Mkn 848 spectral energy density (SED) is lower than 10% (see also Rupke & Veilleux 2013).

X-ray emission has been studied for instance by Iwasawa et al. (2011) and Brightman & Nandra (2011a,b). Chandra data reveal X-ray emission from the two galaxies, and the NW nucleus is much brighter (0.9 dex) than the SE nucleus (Iwasawa et al. 2011). The two systems are instead not resolved in the (deeper) XMM-Newton observations (Brightman & Nandra 2011a), from which more detailed information can be derived for the NW system (assuming that most of the X-ray flux comes from the NW regions; Iwasawa et al. 2011). From the XMM-Newton spectroscopic analysis we derive an observed L_{0.5 − 2 keV} = 1.6 × 10⁴¹ erg s⁻¹; this soft-band emission can be entirely explained by stellar processes (for an SFR = 80 M_⊙ yr⁻¹ we can derive an expected L_{0.5 − 2 keV} ≈ 2 × 10⁴¹ erg s⁻¹; Mineo et al. 2013). Instead, the hard-band shows an excess emission (observed L_{2 − 10 keV} = 3 × 10⁴¹ erg s⁻¹ with respect to ≈10⁴¹ erg s⁻¹ expected from the SFR, Lehmer et al. 2010, although the two values are statistically consistent). The XMM hard spectrum is better reproduced with a heavily obscured AGN, with N_H = 1.4 × 10²⁴ cm² and intrinsic luminosity L_2 − 10 = 3 × 10⁴² erg s⁻¹; the lower confidence limits (c.l.) are of the order of 0.4 dex, while the upper c.l. are unconstrained by the data, which might indicate even higher column densities and intrinsic luminosities. A bolometric luminosity of L_bol = 3 × 10⁴³ erg s⁻¹ can be inferred for the NW SMBH by applying an X-ray bolometric correction of 10 (Lusso et al. 2012). With the currently available archival Chandra and XMM observations, we are unable to derive stringent constraints on the SE X-ray (and bolometric) luminosity and cannot distinguish between dominating AGN and stellar activity in this system. Future observations will allow us to better explore the potentially dual AGN nature of Mkn 848 through hard X-ray emission.

The nuclear optical emission has been studied in detail by Rupke & Veilleux (2013): they observed the NW and SE nuclear regions (white boxes in Fig. 1, right panel), with the Integral Field Unit (IFU) GMOS at the Gemini North telescope. With total exposure times of 2 h (SE) and 2 h 30 min (NW), they detected high-v, blue-shifted ionised (Hα and [N II]λλ6549,6585) and atomic (Na ID) gas associated with conical outflow in the NW system, and the possible presence of ionised outflow in the SE galaxy. Both nuclei were classified as SB-AGN composite systems, and the observed outflows were associated with SF-driven winds. The availability of the deeper MaNGA observations, which have a larger field of view (FOV; dashed line in Fig. 1, right panel) and larger wavelength range than the data presented by Rupke & Veilleux (2013), allowed us to perform a more detailed analysis of the system.

We also note that archival Mkn 848 mm observations could also suggest the presence of molecular outflows. Carbon monoxide CO(3-2) line observations have been obtained by Leech et al. (2010) with the James Clerk Maxwell Telescope. These observations covered the two nuclear regions with a single-pointing and revealed a well-detected, narrow CO profile, with an FWHM ∼ 120 km s⁻¹, and the possible presence of faint extended wings at both negative and positive velocities of up to ≈ ± 250 km s⁻¹. CO(1-0) and CO(2-1) observations have previously also been reported (Sanders et al. 1991; Chini et al. 1992, respectively), but their lower S/N does not permit to confirm or rule out the presence of such high-v wings (see also Yao et al. 2003).

These are some reasons why Mkn 848 can be considered as an ideal laboratory for the study of the SB-QSO galaxy evolution paradigm. In the next sections we present the MaNGA spectroscopic analysis we performed to fully characterise the dynamical and physical conditions in this system.

3. Spectroscopic analysis of the MaNGA data

MaNGA is a new SDSS survey designed to map the composition and kinematic structure of 10⁴ nearby galaxies using IFU spectroscopy. Its wavelength coverage extends from 3600 to 10 400 Å, with a spectral resolution of ∼2000. The diameters of MaNGA data-cubes range from 12″ to 32″, with a spatial sampling of 19–127 fibers. Mkn 848 (MaNGA ID 12-193481) has been targeted with 127 fibers, which almost entirely covers the merging system. The spectroscopic analysis presented in this work takes advantage of the automatic fully reduced MaNGA data release (we refer to Law et al. 2016 for the MaNGA data reduction pipeline).

We performed the data-cube spectral fit following the prescriptions in our previous studies (e.g. Perna et al. 2017a) and adapted them here to trace the kinematic properties of the neutral and ionised ISM components of Mkn 848. Prior to modelling the ISM features, we used pPXF routines (Cappellari & Emsellem 2004; Cappellari 2017) to subtract the stellar contribution from continuum emission and absorption line systems, using the MILES stellar population models (Vazdekis et al. 2010). The separation between stellar and ISM contributions is required in particular to properly measure the kinematics of neutral gas in the ISM. The pPXF fit was performed on binned spaxels using a Voronoi tessellation (Cappellari & Copin 2003) to achieve a minimum S/N >  50 on the continuum in the [5380 ÷ 5770] Å rest-frame wavelength range, modelling the stellar features in the wavelength range 3550–7400 Å. A pure emission and absorption ISM line data-cube was therefore created by subtracting the best-fitting pPXF models.

3.1. Simultaneous modelling of optical ISM features

We used the Levenberg–Markwardt least-squares fitting code MPFITFUN (Markwardt 2009) to simultaneously reproduce the ISM line data-cube with Gaussian profiles, that is, the Hβ, [O III]λλ4959,5007 doublet and He Iλ5876 emission lines, the Na ID absorption line doublet at λλ5891,5897, and the [O I]λ6300, Hα, [N II]λλ6549,6585, and [S II]λλ6718,6732 emission lines.

As mentioned above, atomic and ionised outflowing material has been observed in the NW component of Mkn 848: prominent blue wings have been well detected in the rest-frame optical emission lines and in Na ID absorption features, tracing high-velocity approaching warm and cold material, with ∼ − 900 km s⁻¹ and ∼ − 1200 km s⁻¹, respectively. We therefore modelled these line features as described below.

3.1.1. Ionised lines

We used three (at maximum) sets of Gaussian profiles to reproduce the asymmetric ionised emission line features:

NC. One narrow component (NC) set: 10 Gaussian lines, one for each emission line (Hβ, the [O III] doublet, He I, [O I], Hα, and the [N II] and [S II] doublets) to account for unperturbed systemic emission. The width (FWHM) of this kinematic component was set to be ≲400 km s⁻¹.

OC. One (two) outflow component (OC) set(s): nine Gaussian lines, one for each emission line (without He I, which is usually faint and never well constrained) to account for outflowing ionised gas. No upper limits were fixed for the width of this kinematic component.

We constrained the wavelength separation between emission lines within a given set of Gaussian profiles using the rest-frame vacuum wavelengths. Moreover, the relative flux of the two [N II] and [O III] doublet components was fixed to 2.99 and the [S II] flux ratio was required to be within the range 0.44 <  f(λ6716)/f(λ6731) < 1.42 (Osterbrock & Ferland 2006). All emission line features were fitted simultaneously to significantly reduce the number of free parameters and thus the model degeneracy.

3.1.2. Neutral lines

The sodium absorption lines were fitted simultaneously to the emission features described above. However, a different approach was required to model these lines. Na ID is a resonant system, and both absorption and emission contributions can be observed across the extension of a galaxy (e.g. Rupke & Veilleux 2015; Prochaska et al. 2011). We note that sodium emission is not a common feature: it has been observed at very low level in stacked SDSS spectra (Chen et al. 2010; Concas et al. 2017) and in a few nearby galaxies (see Rupke & Veilleux 2015 and references therein). Unconstrained in the previous observations of Mkn 848, it is clearly detected in the MaNGA data-cube (see Fig. 6).

The modelling of sodium resonant transitions is made all the more difficult because the Na ID system is generally faint, and both systemic and outflowing gas can be associated with emission and absorption contributions (e.g. Prochaska et al. 2011). To reproduce the sodium absorption lines, three (at maximum) sets of Gaussians were used:

NC^±. One narrow component: two Gaussian lines, one for each Na ID transition, to account for (emitting or absorbing) systemic gas. Gaussian amplitudes can vary from negative to positive values to account for absorbing or emitting unperturbed gas. Their kinematic properties (i.e. line widths and centroids) were set to be equal to those of the NC emission lines.

$O{C}_1^{-}$. One outflow component for absorbing gas: two Gaussian lines, one for each emission line of the doublet, to account for outflowing absorbing neutral gas. No upper limits were fixed for the width of this kinematic component, and no constraints were set for the Gaussian centroids (but their relative separation was fixed).

$O{C}_2^{-}$. A second outflow component for absorbing gas: two Gaussian profiles, to account for the presence of very high velocity outflowing gas. Their line widths were fixed to values higher than 700 km s⁻¹; no constraints were set for the Gaussian centroids (but, as before, their relative separation was fixed).

OC⁺. One outflow component for emitting gas: one Gaussian profile, to account for the possible presence of faint wings in emitting Na ID profile. No constraints were set for this Gaussian profile. In our procedure, this component was always fitted together with the sodium OC⁻ Gaussian profiles.

The reduced number of Gaussian profiles in the OC⁺ (with respect to the absorption) is justified by the necessity of minimising the degeneracy in best-fit models; in addition, we note that this broad emission component is generally very faint (see e.g. Fig. 6). Instead, the two Gaussian profiles used for the negative OCs are generally associated with stronger components (Fig. 2).

Fig. 2. Portion of MANGA integrated spectra of Mkn 848, for the NW and SE regions, with superimposed best-fitting components. Top panels: [O III]+Hβ emission lines from the NW (left panel) and SE (right panel) regions. Central panels: He I emission line and Na ID absorption line system. Bottom panels: [N II], Hα and [S II] lines. The different sets of Gaussian components are shown with solid green curves (NC), dashed and dot-dashed blue lines (OC). Stellar continuum emission has been previously modelled with the pPXF routine and subtracted. The dotted lines mark the wavelengths of Hβ and [O III]λλ4959,5007, in the top panels, He I and Na ID doublet lines in the central panels, and [N II]λ6548, Hα, [N II]λ6581 and [S II] doublet lines in the bottom panels. The systemic redshift is defined on the basis of the best-fit solution for the NC in the NW integrated spectrum. The SE systemic lines show a velocity offset of ∼90 km s−1. The best-fits reproduce all but the Na ID NC in the NW integrated spectrum with high significance.

Because the Na ID system consists of two transitions of the same ion, the flux ratio of the f(λ5890) and f(λ5896) lines (hereinafter referred to as the H and K features, respectively) can vary between the optically thick (f(H)/f(K) = 1) and thin (f(H)/f(K) = 2) limits (e.g. Bechtold 2003; Rupke & Veilleux 2015); we therefore constrained the flux ratio within the range 1 <  f(H)/f(K) < 2 for all NC^±, ${\mathrm{OC}}_1^{-}$, and ${\mathrm{OC}}_2^{-}$ Gaussians.

3.1.3. Bayesian information criterion

We performed each spectral fit three times at maximum, with one to three Gaussian sets (e.g. NC, and one or two OC sets for the ionised component) to reproduce complex line profiles where needed. The number of sets used to model the spectra was derived on the basis of the Bayesian information criterion (BIC; Schwarz 1978), which uses differences in χ² that penalise models with more free parameters. In more detail, for an individual model we computed the BIC as χ² + k ln(N), with N the number of data points and k the number of parameters of the model. For an individual spectrum, we first derived ΔBIC_1 − 2 = BIC₁ − BIC₂, where BIC₁ and BIC₂ are the values derived from the models with one and two sets of Gaussians, respectively. When ΔBIC_1 − 2 was lower than a given threshold (ΔBIC_1 − 2 <  30; Mukherjee et al. 1998; Harrison et al. 2016), we favoured the fit with the lower number of components, that is, the one-component fit; otherwise, we performed an additional fit with three components and favoured the two- or three-component models according to the ΔBIC_2 − 3 value, using the threshold described above.

3.2. Non-parametric velocity analysis

The simultaneous line-fitting technique introduced in the previous sections allows us to distinguish between unperturbed gas components (NC), which are associated with emitting gas under the influence of the gravitational potential of the system, and the OC, which most probably is associated with emission of outflowing gas. Therefore, this modelling allows a good separate estimate of the unperturbed and outflowing gas fluxes. Instead, to characterise the kinematic properties of the ISM gas across the FOV best in a homogeneous way and to avoid any dependence on the number of distinct Gaussians used to model the line features in individual spaxels, we used the non-parametric approach (e.g. Zakamska & Greene 2014). Non-parametric velocities were obtained by measuring the velocity v at which a given fraction of the total best-fit line flux is collected using a cumulative function. The zero-point of velocity space was defined by adopting the systemic redshift derived from the NW integrated spectrum (Sect. 4.1). We carried out the following velocity measurements on the best fit of [O III]λ5007 ([O III] hereinafter) and Hα line profiles for the ionised gas and of the Na ID to trace the neutral ISM component:

• W80: The line width comprising 80% of the flux, which is defined as the difference between the velocity at 90% (V₉₀) and 10% (V₁₀) of the cumulative flux.

• V₀₂: The maximum blue-shift velocity parameters, which is defined as the velocity at 2% of the cumulative flux.

• V₉₈: The maximum red-shift velocity parameters, which is defined as the velocity at 98% of the cumulative flux.

• V_p: The velocity associated with the emission line peak.

These velocities are usually introduced to characterise asymmetric line profiles. In particular, W80 is closely related to the velocity dispersion for a Gaussian velocity profile, and is therefore indicative of the line widths, while maximum velocities are mostly used to constrain outflow kinematics when prominent line wings are found (e.g. Cresci et al. 2015a; Perna et al. 2017a; see also e.g. Harrison et al. 2018). We used the 2% and 98% velocities to facilitate comparison with Rupke & Veilleux (2013), although the authors defined the maximum blue-shifted velocity (V₉₈ in their work) in a slightly different way and derived the outflow velocity from the broad (outflow) component alone. As a result, their estimates may be higher than our values (which are defined from the total emission line profile). Finally, V_p, which generally traces the velocity of the brightest narrow emission line components, can been used as a kinematic tracer for rotating-disc emission in the host galaxy (e.g. Harrison et al. 2014; Brusa et al. 2016; Perna et al. 2018a).

We note that perturbed kinematics that are due to gravitational interactions of the two galaxies could also be present in the system, especially along the tidal features. We used our multi-component models and the non-parametric techniques to distinguish between gravitational interaction and rotational and outflow signatures.

4. Cool and warm gas kinematics

In this section we present the results we obtained from the fitting prescription and the non-parametric approach exposed above. We show integrated spectra with high S/N of the two merging galaxies and the spatially resolved kinematic and physical properties of Mkn 848.

4.1. Integrated spectra

Before analysing the spatially resolved information of the ISM in detail, we extracted the integrated NW and SE nuclear spectra from the MaNGA data-cube, considering the spaxels associated with S/N >  45 in the stellar continuum at ∼5550 Å. The two spectra are shown in Fig. 2. The NW and SE galaxies are at z = 0.0404 and 0.0407, respectively. In the following analysis, we derived all velocities using the systemic of the NC in the NW integrated spectrum as a reference unless stated otherwise.

Asymmetric profiles with extended blue and red wings are observed in all emission lines. In particular, OC Gaussians are blue-shifted with respect to the NC in the NW region; in the SE spectrum, all the OCs are instead red-shifted with respect to the narrow Gaussians. The red-shifted emission in the SE nucleus was not detected by Rupke & Veilleux (2013), probably because of the low fluxes associated with this nucleus (f_SE : f_NW ∼ 1 : 10 in all the emission line features; Fig. 2). We also note that all the profile peaks are associated with the NC Gaussians, and that outflowing gas therefore does not dominate in the emission of ionised gas.

The sodium absorption systems display very complex profiles. In the NW spectrum, the NC are faint and poorly constrained (but their kinematic properties are fixed to those of the other emission lines), and the bulk of the absorption is associated with outflowing gas. The SE spectrum displays similar complexities, with prominent outflow components. In this case, the neutral outflow in the SE system was not previously detected with the GMOS instrument either (Rupke & Veilleux 2013).

Taking advantage from the spatially resolved analysis presented in the next sections, we also show in Fig. 3 the outflow-integrated spectra, which were obtained considering the only spaxels associated with very high-velocity in the neutral and/or ionised phases. These spectra highlight the high-velocity blue and red wings in the ionised gas line profiles, and more importantly, the complexities in the emission and absorption line profiles. In particular, the [O III] and Na ID outflow tracers were fitted with two or three Gaussian components. The very broad and smoothed profiles do not allow us to uniquely identify the individual kinematic components associated with the outflows for either the ionised or atomic phases.

Fig. 3. NW (left panel) and SE (right panel) outflow-integrated spectra extracted from the spaxels associated with high-velocity gas in the neutral and/or ionised phases (ΔV02([O III]) <  −550 km s−1 and/or ΔV02(Na ID) <  − 900 km s−1, for the NW system; ΔV98([O III]) >  450 km s−1 and/or ΔV02(Na ID) <  − 750 km s−1, for the SE galaxy; see Sect. 4.2). The shell-like region spaxels are not included in these integrated spectra. The stellar continuum emission has been previously modelled and subtracted from each spaxel. The red curves overplotted on the spectra represent the best-fit solutions that reproduce the line profiles; NC and OC Gaussians are also shown as in Fig. 2. The outflow-integrated spectra extracted from the original MaNGA datacube (hence without pPXF best-fit model subtraction) and the integrated pPXF best-fit spectra are also shown in the figure with solid and dotted grey curves, respectively. For the two grey curves we applied a downward shift for visualisation purposes.

In Table 1 we report the spectroscopic results obtained from modelling the NW and SE integrated spectra. In the next sections, we introduce different line ratio diagnostics to constrain the ISM conditions (e.g. dust extinction and ionisation; see Table 1). In advance of this, we note that the flux ratios f([O III])/f(Hβ) and f([N II])/f(Hα) that are normally used to separate purely SF galaxies from galaxies containing AGN (Baldwin et al. 1981) suggest a composite SB-AGN ionisation in the two systems. The outflowing gas is instead associated with higher flux ratios that are compatible with AGN ionisation (see Sect. 8 for details).

4.2. Spatially resolved spectroscopy

To derive spatially resolved kinematic and physical properties of ionised and atomic gas, we analysed the spectra of individual regions in the FOV. Before proceeding with the fit, we smoothed the data-cube using a Gaussian kernel with an FWHM of three pixels (1.5″), taking into account the seeing conditions at the time of MaNGA observations (between 1.3″ and 1.9″). A Voronoi tessellation was also derived to achieve a minimum S/N = 10 of the [O III]λ5007 line for each spaxel. We then applied BIC selection in each Voronoi bin to determine where a multiple-component fit was required to statistically improve the best-fit model. This choice allowed us to use the more degenerate multiple-component fits only where they were really needed. In order to mitigate the degeneracies in our fit as well as possible, we also considered an additional constraint for the relative fluxes of emission line features, requiring all the OC amplitudes to be smaller than those of NC (per individual line). This constraint was adopted because of the best-fit results obtained from the nuclear and outflow-integrated spectra; a visual inspection of individual fits and related residuals confirmed the advantages of using this.

Figure 4 (left) shows a map reporting the number of sets of Gaussians we used to model the spectra in individual spaxels. As expected, most of the two- and three-set models were selected for the spaxels with the highest S/N, that is, in the regions close to the two nuclei¹. The figure also displays the [O III]λ5007, Hβ, and Hα flux maps obtained by integrating the best-fit total emission line profiles. Some off-nuclear line emission is detected in the NW tidal tail and in the lower part of the figure, where we probably observe the overlapping of the two tidal tails that we identified in the HST images (T1 and T2 in Fig. 1, third panel).

Fig. 4. Top-left panel: number of sets of Gaussians we used to model individual spaxel spectra. Top-right and bottom panels: Hβ, [O III], and Hα integrated fluxes obtained from simultaneous best fits (we show only the spaxels with an S/N >  3). Solid contours refer to S/N levels of 10 and 40 in the stellar continuum map (at ∼5550 Å). Black crosses indicate the nuclear regions.

Before introducing the velocity characterisation of the merging system, we show in Fig. 5 a wavelength-collapsed image of the Mkn 848 high-v [O III]λ5007 gas (absolute velocities > 250 km s⁻¹ in the NW nucleus reference frame). We also display a number of spectra extracted from 3 × 3 spaxel regions. The spectra associated with the more nuclear regions display profiles that are more highly asymmetric, with prominent blue- and red-line wings, consistent with what we found in the integrated spectra (Sect. 4.1). The extracted spectra also highlight the variations in line peak wavelength across the FOV, and the possible presence of a double-peaked emission line at the overlapping of the T1 and T2 tidal tails (panel f). Broad profiles are found in the eastern quadrant above the SE nucleus, while narrower features are found in the spectra extracted from the T1 tail and the southern elongation.

Fig. 5. High-velocity [O III]λ5007 emission line flux across the FOV obtained by integrating for each spaxel the spectra in the regions 4890 − 5004 Å and 5014 − 5030 Å, which correspond to absolute velocities > 250 km s−1 with respect to the reference frame of the NW nucleus. The map was obtained after subtracting the stellar contribution (using pPXF best-fit spaxel spectra). Solid contours refer to the S/N = 10 and 40 levels in the stellar continuum map. Different spectra, extracted from the regions labelled in the flux maps, are also shown in the insets, with normalised fluxes; blue and red spectra in inset (b) show the most asymmetric line profiles and are extracted from the 3 × 3 pixel regions shown in the central panel with blue and red squares, respectively. Vertical dashed lines indicate the (NW nucleus systemic) wavelength of the [O III]λ5007 line; dotted lines refer to ±250 km s−1 with respect to the NW nucleus systemic velocity.

Figure 6 shows the Na ID absorption and emission line fluxes across the FOV, obtained by integrating the spectral regions 5860 − 5910 Å and 5890 − 5920 Å, respectively, after ionised ISM and stellar contributions were subtracted. The maps show that the emission line contribution is small and is generally found at larger distances from the nuclear regions. We also display different spectra extracted from distinct regions, in which Na ID is dominated by absorbing or emitting contributions.

Fig. 6. Central panels: Na ID absorption (Abs.) and emission (Em.) line fluxes across the FOV obtained by integrating for each spaxel the spectra in the regions 5860 − 5910 Å and 5890 − 5920 Å, respectively. Negative and positive fluxes were obtained after subtracting the stellar contribution (pPXF best fit); for the absorption fluxes, we also removed the He I emission line (simultaneous best-fit results). Insets: different spectra extracted from the regions labelled in the flux maps, with normalised fluxes: absorption- (emission-) dominated Na ID line profiles are shown in the right (left) panels. Red curves represent the Na ID stellar contribution fitted with pPXF, rescaled for visual purposes. Vertical dashed lines indicate the (NW nucleus systemic) wavelengths of the He I and Na ID doublet lines.

In the next sections we show that simultaneous fits successfully reproduce all the kinematic and flux variations described above for the ionised and atomic components. We anticipate, however, that the data quality and the complex broad profiles associated with the atomic and ionised gas tracers do not allow us to distinguish between the different outflow kinematic components that were used to reproduce the line profiles.

The pPXF analysis allowed us to trace the stellar velocities almost throughout the MaNGA FOV (Fig. 7, left). The stellar velocity map displays a strong gradient, from higher to lower velocities extending from SE to NW, taking as a reference the systemic of the NW nucleus; an inversion of the gradient is instead observed along the T1 tidal tail, which probably extends to the position of the SE nucleus (see the HST F160W image in Fig. 1). The extension of the T1 tail near the SE nucleus seems to have a stellar velocity that is similar to the velocities of the SE nucleus (see also Fig. 5).

Fig. 7. Left panel: stellar velocity map derived from pPXF best-fit models. Central panel: [O III] non-parametric peak velocity Vp estimate (an identical map is obtained for the Hα line). Right panel: difference between [O III] peak velocity (Vp) and stellar velocity (Vstar); the most significant discrepancy between the two kinematic tracers is found in regions A and B in the central panel. All velocities are derived taking the narrow-component velocity of the NW -integrated spectrum as a reference.

Similar velocity maps are obtained considering the V_p of the [O III]λ5007 and Hα lines (Fig. 7, central panel), suggesting that the velocity peak of ionised gas emission lines is a good tracer of the stellar motions, as a result of rotation in the individual galaxies and of the gravitational interaction (see e.g. the tidal tail T1). Because of the adopted constraints (Sect. 4.1), the velocity peak is always associated with the NC emission; therefore, we can reasonably assess that the OC does trace outflowing gas.

Significant differences are found between V_p and V_star in the NE and SW regions around the SE nucleus, however (areas A and B in the central panel). This discrepancy, highlighted in Fig. 7 (right), is due to the significant variations in the V_p velocities: in region A, the velocity shift can be explained by the interaction of tidal tails gas (see Fig. 5), while in region B, different kinematics may play a role (e.g. collimated outflow). In the following, we therefore use the stellar velocity map to subtract the gravitational motions and unveil the outflow kinematics.

The V_p velocity maps in the vicinity of the nuclei are similar to those reported by Rupke & Veilleux (2013), who found a rotation pattern in the SE region and strong tidal motions and a near face-on orientation in the NW regions. In the next sections we discuss the kinematic properties of ionised and atomic gas by searching for deviations from the general motions traced by the stellar component and the bulk of ionised gas. In particular, we refer to the non-parametric maximum velocities (V₀₂ and V₉₈) of [O III] and Na ID features to trace the ejected gas kinematics.

4.2.1. [O III] as a tracer of ionised outflow

Because of the strong velocity trend in the merging system, the maximum velocity maps mostly reproduce the same behaviour as the stars and the bulk of ionised gas (Fig. 7). In order to unveil the outflows in the system, we therefore plot in Fig. 8 the difference between the non-parametric maximum velocities and V_star.

Fig. 8. Maps across the FOV of the maximum blue- (left panel) and red-shifted (central panel) velocities for the [O III] best-fit line profile. In particular, the two panels show the differences between maximum velocities (V02 and V98) and stellar velocity (Vstar), as labelled, to reveal outflow signatures. In the V98 − Vstar map the shell region (S) associated with both receding and approaching high-v gas is shown. The right panel shows the [O III] velocity width W80.

In the V₀₂ − V_star map, velocities as high as −550 ÷ −650 km s⁻¹ are observed in the two nuclear regions, with projected extensions of ∼5 kpc. The red-shifted maximum velocities (V₉₈ − V_star) instead display fast-receding gas in the proximity of the SE nucleus, with velocities of the order of ∼450 ÷ 550 km s⁻¹. Signatures of disturbed velocities are also found in the proximity of the region associated with the shell (S), with approaching as well as receding fast gas. In general, the distribution of high-v gas suggests nuclear outflows in the two systems (see also Sect. 4.2.3).

The Hα velocity map confirms the presence and extension of fast gas in the proximity of the two nuclei. We note, however, that although the Hα line is stronger than that of doubly ionised oxygen, it might suffer from blending with [N II] emission lines. This means that the kinematics associated with Balmer emission might be affected by larger uncertainties.

4.2.2. Na ID as a tracer of neutral outflow

The Na ID negative maximum velocities were derived using the wavelength of the H transition of the NW integrated spectrum as a zero-point. Figure 9 shows the maximum velocity maps derived for the sodium absorption line. The negative V₀₂ − V_star tracer (ΔV₀₂ hereinafter for simplicity) highlights strong neutral outflows in the NW nucleus (with ΔV₀₂ ∼ −1100 km s⁻¹), and of slower winds with ΔV₀₂ ∼ −750 km s⁻¹ in the other nucleus.

Fig. 9. Sodium-line velocity maps: the maximum velocity ΔV02 (left panel), highlighting the receding-absorption components, and ΔV98 (central panel), highlighting the approaching-emitting components. The right panel shows the absorption component W80 velocity map.

Because of the geometrical configuration that is required to observe absorption features, we do not expect to observe the neutral counterpart of the red-shifted ionised outflows in absorption lines. On the other hand, red-shifted emission from the same transitions, as seen in the previous section, is detected throughout the FOV: this positive flux might be associated with the receding part of the outflow that is not visible in absorption, in the absence of significant absorption on the line of sight (LOS; Prochaska et al. 2011).

The central panel in Fig. 9 shows the V₉₈ − V_star (ΔV₉₈ hereinafter) velocities derived for the sodium emission features. The highest velocities (∼900 km s⁻¹) are associated with the NW circum-nuclear regions; we note, however, that the low S/N (see Fig. 6) does not allow us to derive robust kinematics for this gas component. Moreover, as described in the previous section, modelling absorption and emission contributions simultaneously leads to significant degeneracy². The simultaneous fit of all the emission and absorption lines allows us to mitigate this degeneracy: the sodium emission line of low-velocity gas is constrained to the same systemic as the ionised gas. We also note that the ΔV₉₈ map qualitatively confirms that our fit models are able to mitigate the degeneracies in the Na ID, showing that the regions associated with very high-velocities are the same where ionised and neutral (absorbing) outflow material is also detected.

MaNGA data allowed us to clearly detect high-v sodium emission in Mkn 848 (e.g. Fig. 6). Resonant line emission from cool receding outflow gas has been observed in stacked spectra and in a handful of individual spectra (e.g. Chen et al. 2010; Rubin et al. 2011; see also e.g. Martin et al. 2013 for Mg II resonant transition). To our knowledge, a spatially resolved resonant emission line from galactic wind has only been observed in two other systems: NGC 1808 (Phillips1993) and F05189-2524 (Rupke & Veilleux 2015).

4.2.3. Multi-phase outflow configuration

Figures 8 and 9 display disturbed kinematics in the proximity of the NW and SE nuclei, as well as in the shell observed between the two systems. The disturbed kinematics involving a multi-phase medium is also highlighted by the [O III] and Na ID (absorption) W80 velocities that are reported in the same figures.

Approaching outflowing ionised gas with velocities as high as ∼ − 650 km s⁻¹ is observed in the NW regions. Neutral ejected material with maximum velocities of ∼ − 1100 km s⁻¹ is also detected in the same regions. These values, in particular those related to the neutral component, are greater than the velocities generally found for SB-driven winds in local star-forming galaxies and ULIRGs (a few hundreds km s⁻¹; e.g. Cazzoli et al. 2016; Chen et al. 2010; Rupke & Veilleux 2013), although this high-v might be explained by smaller projection effects with respect to those of other SB-driven outflows. Instead, they are similar to the velocities found in AGN hosts (e.g. F06206-6315, Cazzoli et al. 2016; Mrk 273 and Mrk 231, Rupke & Veilleux 2013). It is therefore possible that this material is ejected by AGN-driven winds, or at least SMBH activity might play a key role in this process. The presence of receding neutral high-v gas is also detected in these nuclear regions (Fig. 9, central panel).

In the SE system we observe outflowing ionised emission associated with both approaching and receding material, with velocities of ±500 km s⁻¹. The spatial extent of this outflow and the near edge-on orientation of the SE galaxy may suggest that the observed outflow fills a wide biconical region that probably encloses the high-v gas close to the shell-like region. Approaching and receding neutral outflows, reaching velocities of ΔV_max ≈ ±700 km s⁻¹, are instead only revealed in the more nuclear regions because of observational limitations: the detection of absorbing Na ID is related to strong background emission.

Our spatially resolved spectroscopic analysis allows us to detect high-v ionised gas out to a projected distance of $R_{\mathrm{out}}^{\mathrm{NW}}([\mathrm{O}\mathrm{III}])\approx 4$ kpc for the outflow associated with the NW nucleus, and of $R_{\mathrm{out}}^{\mathrm{SE}}([\mathrm{O}\mathrm{III}])\approx$ 6 kpc for the one related to the SE system (see Fig. 8). The observational limitations affecting the Na ID analysis allow us to derive lower limits on the real neutral gas widening, with $R_{\mathrm{out}}^{\mathrm{NW}}(\mathrm{Na}\mathrm{ID})\gtrsim 4$ kpc and $R_{\mathrm{out}}^{\mathrm{SE}}(\mathrm{Na}\mathrm{ID})\gtrsim 2$ kpc (Fig. 9).

In Sect. 9 we use these distances to derive outflow energetics assuming that the outflows extend out to these radii from the central regions. We note that the near edge-on orientation of the SE galaxy suggests that the de-projected distances are not so different from the estimates reported before. On the other hand, the almost face-on configuration of the NW galaxy does not allow us to constrain the related outflow extension well, and large uncertainties might therefore be associated with these estimates.

While the geometrical configuration of the SE galaxy and the detection of multi-phase approaching and receding gas strongly suggest the presence of a biconical outflow whose axis is perpendicular to the galaxy major axis, the NW system could be more complex. Its almost face-on orientation would be compatible with a biconical outflow whose red-shifted gas component is mostly obscured by the host galaxy. Nevertheless, we traced faint high-v Na ID gas emission, which is associated with receding gas. This cool material might suggest that the outflow in the NW galaxy has a wide aperture, such that the far parts of the cone that intercept the plane of the sky, are associated with positive velocities (e.g. Perna et al. 2015a; Venturi et al. 2018b). Alternatively, this system might be associated with a more complex configuration between the outflow axis, the galaxy, and the plane of the sky (see e.g. Fischer et al. 2013). Data with a higher spatial resolution and S/N are required to distinguish between the different scenarios; for the purposes of this work, we assume that the two merging galaxies are associated with biconical outflows extending out to the distances R_out reported above (but see Rupke & Veilleux 2013 for an alternative configuration of the NW outflow).

5. Dust absorption

Another key ingredient for characterising the outflow properties is the obscuration related to the ejected material. We used the Balmer decrement ratios f(Hα)/f(Hβ) to estimate the reddening across the FOV. A ratio of 3.1 (Case B, Osterbrock & Ferland 2006) is mostly used to determine the amount of extinction for low-density gas, such as that of the narrow-line regions (NLR). Any deviation to higher flux ratios is interpreted as due to dust extinction, which obscures emission lines at lower wavelengths.

Figure 10 shows the flux ratios separately for the total Balmer profiles and for NC and OC. The higher flux ratios, corresponding to the higher extinctions, are related to the nuclear positions and to the dust lane in the SE galaxy that is also visible in the three-colour HST image (Fig. 1). The OC components are associated with even higher dust extinctions in the regions associated with the highest outflow velocities (Fig. 8). This is consistent with the measurements reported in the literature for AGN-driven outflows (e.g. Holt et al. 2011; Mingozzi et al. 2019; Perna et al. 2015a; Villar-Martín et al. 2014). However, such high dust extinctions might be the result of an incorrect subtraction of the stellar component, which would mostly affect the Hβ emission line wings (see e.g. Perna et al. 2017a). In Fig. 3 we show the MaNGA spectra and the best-fit pPXF spectra integrated over the outflow regions in the NW and SE galaxies. These integrated spectra show that the pPXF best fit does on average not underestimate the fluxes below the Hβ wings (the best fit and the integrated spectra match almost perfectly, except for a negligible discrepancy in the red wings). We can therefore reasonably exclude that the high dust extinctions associated with the OC are due to an incorrect subtraction of the stellar component.

Fig. 10. Balmer decrement f(Hα)/f(Hβ) ratio across the FOV for the total line profiles (left panel) and for the NC (centre) and OC (right panel) components separately. Higher flux ratios correspond to higher dust extinction.

To explain these high extinctions, two possible explanations have been proposed: (i) the outflows might be obscured by a dusty natal cocoon, or (ii) ionised and dusty phases are mixed in the outflowing gas. Flux ratios above the 3.1 value are measured thorughout the entire FOV, which confirms the dust-obscured nature of the interacting system, but the highest flux ratios are spatially correlated with fast outflowing material. This argument, together with the fact that we observe both ionised and atomic outflows, points toward the scenario in which warm and dusty outflowing phases are mixed.

For the Case B ratio of 3.1 and the SMC dust-reddening law, the measured Balmer decrement ratios translate into a V-band extinction A_V in the range between 1.3 and 3.1. These values are consistent with those derived from the integrated nuclear and outflow spectra presented in Sect. 4.1 (Table 1), although the integrated measurements cannot account for the variations observed in Fig. 10.

6. Cool gas absorption and hydrogen column densities

In this section we derive the column density (N_H) of the wind, which is required to estimate how much neutral material is brought out in the outflow. The column density can be derived using the so-called doublet-ratio approach, which relates N_H with the sodium equivalent width via the sodium curve of growth (e.g. Cazzoli et al. 2016). This approach assumes that the absorbing neutral material has a low optical depth (τ ≲ 1) and a uniform and total coverage of the continuum source, with a covering factor C_f = 1. More plausible scenarios instead consider distributions of τ for individual kinematic components that also extend from low to high optical depths and C_f ≠ 1 (e.g. Rupke et al. 2005a and references therein), as well as the possibility that sodium atoms at different velocities might be exposed to a different background continuum source (because of overlapping clouds on the LOS; see Fig. 6 in Rupke et al. 2005a).

The very complex Na ID profiles and the low spectral resolution and quality of MaNGA data do not allow us to derive reliable spatially resolved column densities, especially because of the velocity overlap between the different kinematic components in the sodium profiles. Moreover, the multi-Gaussian fit performed in our analysis to derive the sodium kinematics instead assumes that all the kinematic components are co-spatial along the LOS and are exposed to the same background continuum source.

Because of these limitations, we decided to derive the average column densities associated with the neutral outflows in the NW and SE systems from the integrated spectra extracted from the spaxels associated with high-velocity gas in the neutral phase (see Sect. 4.1). These spectra, which are characterised by higher S/N, allowed us to follow a more rigorous approach. We fitted the observed intensity profiles of the sodium lines with a model parameterised in the optical depth space (e.g. Rupke et al. 2002, 2005a),

$$\begin{array}{cc}\hfill I(\lambda )& =1-C_f\times [1-\exp (-{\tau }_0{\mathrm{e}}^{-{(\lambda -{\lambda }_K)}^2/{({\lambda }_{K}b/c)}^2}\hfill \\ \hfill & -2{\tau }_0{\mathrm{e}}^{-{(\lambda -{\lambda }_{\mathrm{H}})}^2/{({\lambda }_{\mathrm{H}}b/c)}^2})],\hfill \end{array}$$(1)

where I(λ)=f(λ)/f₀(λ) is the continuum-normalised spectrum, C_f is the covering factor, τ₀ is the optical depth at the line centre λ_H (5891 Å), b is the Doppler parameter ($b=FWHM/[2\sqrt{ln(2)}]$), and c is the velocity of light. This model assumes that the velocity distribution of absorbing atoms is Maxwellian and that C_f is independent of velocity. Following Rupke et al. (2005a), we assumed the case of partially overlapping atoms on the LOS, so that the total sodium profile can be reproduced by multiple components and $I(\lambda )={\mathrm{\Pi }}_{i=1}^n{I}_i(\lambda )$, where I_i(λ) is the ith component (as given in Eq. (1)) used to model the sodium features (see also Sect. 3.1 in Rupke et al. 2002).

The best-fit models shown in Fig. 11 were obtained assuming one outflow component for the SE system and two outflow components for the NW galaxy, constraining their kinematics (FWHM and line centres) to the values previously obtained from the Gaussian fits (Fig. 3). The existence of a multi-component outflow in the NW system is suggested by negative fluxes at ≈5870, next to the blue wing of the He I emission line. This contribution is also clearly visible in Fig. 6, panel b, and in the SDSS legacy survey spectrum (Perna et al. 2017a). We associate it with the Na ID absorption that also clearly affects the red wing of the He I line (see Figs. 2 and 3). The feature at ≈5870 Å suggests very high-v gas (V_max ≈ −1200 km s⁻¹), which is also confirmed by Rupke & Veilleux (2013, see their Table 3). To obtain the best-fit models shown in Fig. 11, we also constrained the covering factors and optical depths within the limits 0 ≤ τ₀ ≤ 10³ and 0 ≤ C_f ≤ 1.

Fig. 11. Portion of the outflow-integrated spectra of Mkn 848 around the He I and Na ID lines for the NW (left panel) and SE (right panel) regions. The spectra are extracted from the spaxels associated with high-velocity gas in the neutral phase, with ΔV02 <  −900 km s−1, for the NW system, and ΔV02 <  −750 km s−1 for the SE system. The two spectra have a normalised continuum, and the superimposed best-fitting Na ID components in absorption are obtained following Rupke et al. (2005a), assuming a geometrical configuration with partially overlapping sodium atoms (Sect. 6). The He I emission line is instead fitted with a Gaussian profile.

We obtained for the outflow components covering factors of the order of 0.2 − 0.5 and low optical depths (τ₀ ≈ 0.4 − 0.8). The column densities were computed from the best-fit values of τ₀ and b using Eq. (7) in Rupke et al. (2002).

Finally, following the prescriptions reported in Rupke et al. (2005b), we derived the sodium abundance taking into account the ionisation fraction (y = 0.9) and the depletion of Na atoms onto dust and the abundance as

$$\begin{array}{c}\hfill log(N_{\mathrm{Na}}/N_{\mathrm{H}})=log(1-y)+A+B,\end{array}$$(2)

where A = log(N_Na/N_H)_gal is the Na abundance in the galaxy and B = log(N_Na/N_H)−log(N_Na/N_H)_gal = −0.95 is the depletion (the canonical Galactic value, Savage & Sembach 1996). Using the galaxy abundance derived from the Eq. (12) in Rupke et al. 2005b, we derived log(N_H/cm⁻²) = 20 ± 0.1 for the NW and SE outflows³. These values are used to derive the atomic outflow energetics in Sect. 9.

7. Plasma properties

Plasma properties represent further key ingredients to derive outflow energetics. Electron density (N_e) and temperature (T_e) within AGN-driven outflows can be derived using diagnostic line ratios involving forbidden line transitions (Osterbrock & Ferland 2006). At optical wavelengths, the [S II]λλ6716,6731 and [O III] flux ratios (involving [O III]λλ4959,5007 and [O III]λ4363) have been used in the literature to derive the electron density and temperature for a limited number of sources (Perna et al. 2017b and references therein). These estimates are quite uncertain, however: doublet ratio estimates for individual kinematic components are strongly affected by the degeneracies in the fit because the involved emission lines are faint and blended (see e.g. the detailed discussion in Rose et al. 2018). Moreover, the [S II] line ratio is sensitive to relatively low densities (10² ≲ N_e/cm⁻³ ≲ 10^3.5), while ionised outflowing gas has been found also in higher N_e regions (e.g. Lanzuisi et al. 2015; Villar-Martín et al. 2015). Furthermore, the [O III] and [S II] species, being associated with different ionisation potentials, could arise from different emitting regions. However, while different line diagnostics have been proposed in the literature (e.g. involving trans-auroral emission lines; e.g. Holt et al. 2011; Rose et al. 2018), the [S II] line ratios still represent the more feasible diagnostic to derive the outflow electron density: the trans-auroral diagnostics are very faint and can only be detected in local targets; moreover, the use of these diagnostics requires a wide spectroscopic wavelength coverage.

In Fig. 12 we show the [S II] ratio map, derived by computing the line ratio in each spaxel where the sulfur lines are detected with an S/N >  3 for the total line profiles and for the NC and OC separately. The NC map shows that the lower ratios, which are associated with higher electron densities, are found in the proximity of the NW nucleus, with N_e ∼ 150 ÷ 300 cm⁻³ (from Eq. (3) in Perna et al. 2017b, assuming T_e = 10⁴ K). The SE nuclear regions are instead associated with lower electron densities (≲100 cm⁻³).

Fig. 12. [S II] flux ratio maps derived for the total line profiles (left panel), and distinguishing between systemic (centre) and outflowing gas (right panel). Higher flux ratios correspond to lower electron densities.

The outflow line ratios are lower than those of NC, especially in the regions that are associated with the high-v gas in the NW nucleus. They correspond to electron densities as high as N_e ≈ 1400 cm⁻³. The SE nucleus is instead associated with high line ratios that correspond to N_e <  100 cm⁻³.

Because the [O III]λ4363 feature is so faint, we were unable to derive an electron temperature map from the [O III] ratio f(λ5007 + λ4959)/f(λ4363); we therefore used the outflow-integrated spectra obtained in Sect. 4.1 to derive average electron temperatures associated with the outflowing gas. In Table 2 we report the T_e and N_e obtained from the integrated spectra; we note that the electron density values are consistent with the individual spaxel estimates (Fig. 12).

The wide wavelength coverage of MaNGA observations allowed us to derive further estimates of the plasma properties from the NW- and SE-integrated spectra, in which other faint temperature- and density-sensitive lines can be detected. In particular, we used the [S II] and [N II] temperature-sensitive line ratios, the [O II] density-sensitive ratio (Osterbrock & Ferland 2006), and the trans-auroral line diagnostic to derive N_e (see Table 2 for details). We modelled these additional emission lines using a multi-component simultaneous fit; for the features whose S/N was high enough to distinguish between different kinematic components, we derived temperature- and density-sensitive flux ratios for both unperturbed and outflowing ionised gas. All the line ratio diagnostics and the T_e and N_e measurements are reported in Table 2.

We found that regardless of the small discrepancies between the values obtained with different tracers, systemic gas is in general associated with low electron densities (up to a few 100 cm⁻³), while outflowing emission is associated with similar (in the SE regions) or higher densities (in the NW system; $N_{\mathrm{e}}^{\mathrm{NW}}(\mathrm{OC})\approx 1500$ cm⁻³). These results are consistent with previous works, in which both stacked spectra (Perna et al. 2017b) and high S/N spectra of individual targets (e.g. Villar-Martín et al. 2015; Kakkad et al. 2018; Mingozzi et al. 2019; Rose et al. 2018) have been analysed.

Lower electron density values were instead derived from the trans-auroral line diagnostics, for which we were unable to distinguish between systemic and outflow components. The reason probably is that the systemic gas, associated with lower N_e, dominates the total line ratio fluxes. We also note that the oxygen-line doublet, being very close in wavelength, does not allow us to derive reliable measurements for N_e.

Electron temperatures of the order of 2 × 10⁴ K were derived from the [O III] line ratios for the systemic and outflowing gas, consistently with the results presented by Perna et al. (2017b), for the NW and SE systems (Table 2). Lower temperatures were instead measured from the sulfur and nitrogen diagnostics (T_e ≲ 10⁴ K). These discrepancies are typical of photoionised gas and can be explained considering the different excitation potentials of the different transitions: [N II] and [S II] are stronger in the outer parts of the photoionised regions, where the temperature is lower (Osterbrock & Ferland 2006; see Curti et al. 2017). In Sect. 9 we use the [S II]- and [O III]-based outflow plasma properties we derived from the NW and SE outflow-integrated spectra to constrain the energetics of the ejected material in the two systems.

8. Ionisation conditions

In this section we investigate the dominant ionisation source for the emitting gas across the MaNGA FOV using the BPT (Baldwin et al. 1981) diagram. Figure 13 (top panels) shows the [O III]λ5007/Hβ vs. [N II]λ6584/Hα flux ratios, derived by integrating the line flux over the entire profiles (left panel) and the outflowing gas components (right), for only those spaxels in which all the emission lines (or the OC) are detected with an S/N >  3⁴. The curves drawn in the diagram correspond to the maximum SB curve (Kewley et al. 2001) and the empirical relation (Kauffman et al. 2003) used to separate purely SF galaxies from composite AGN-H II galaxies and AGN-dominated systems (e.g. Kewley et al. 2006). The points in the diagnostic diagrams are colour-coded based on the different subset selected in the BPT; the same colours are also reported in the map of Mkn 848 (insets in the figure) to distinguish the different ionisation sources in the FOV.

Fig. 13. Top panels: standard diagnostic [O III]/Hβ vs. [N II]/Hα diagrams. The lines drawn in the figures display the curves used to separate H II-like (below the dashed line) from AGN-ionised regions (above the solid line); the locus in between the two curves is associated with composite SF-AGN regions. The solid line is adapted from Kewley et al. 2001; the dashed curve is taken from Kauffman et al. (2003). In the left panel, flux ratios are obtained considering the total line profiles (i.e. NC+OC fluxes); in the right panel, only the OC fluxes have been used. Small colour-coded dots refer to single spaxel measurements associated with the spatial regions shown in the insets (axis scales are the same as in all previous maps). Black spaxels are associated with flux ratios below the curve from Kauffman et al. (2003); purple to red instead show the increasing distance from the same curve. Central panels: standard diagnostic [O III]/Hβ vs. [S II]/Hα diagrams. The lines separates H II-like ionised regions (below the solid line) from AGN-ionised gas (above the dotted and dashed lines) and LINERs (in the right part of the diagrams). The solid and dashed lines are taken from Kewley et al. (2001) and Kewley et al. (2006), respectively. Colour-coded flux ratio dots highlight the separation between AGN, LINERs, and H II regions, and refer to the spatial regions shown in the insets. Bottom panels: standard diagnostic [O III]/Hβ vs. [O I]/Hα diagrams. The colours of the flux ratio dots and separation lines are the same as those in the central panels. In all figures, the largest blue and red dots refer to the flux ratios obtained from the NW and SE outflow-integrated spectra, respectively; only the flux ratios associated with Hβ and [N II] (or [S II], in the central panels; or [O I] in the bottom panels) lines with an S/N >  3 in are shown in the figures. Representative 1σ errors are shown in the bottom corners of the plots.

The BPT diagram derived from total line profiles shows that the H II-like ionisation is found along the T1 tidal tail, and that the emission in the NW system is mostly associated with composite AGN-H II flux ratios. Both the strong star formation activity and the (possible) presence of a buried AGN (see Sect. 2) in the nuclear region may therefore be responsible for the observed emission line ratios in the NW galaxy. As mentioned in Sect. 4.2.3, the observed outflow velocities can be used as a discriminant of SB- vs. AGN-driven winds, pointing to the latter mechanism. However, it is possible that SB contributes to the outflow. Other diagnostics are therefore required to constrain the nature of the observed atomic and ionised outflows.

The ionisation conditions in the SE system are instead more heterogenous: AGN-H II line ratios are observed only along the dust lane, while the directions perpendicular to the galaxy major axis are associated with AGN flux ratios. The presence of high-v ionised and atomic gas in the same direction (see Figs. 8 and 9) allows us to argue that the disturbed kinematics in the SE system are probably due to an AGN-driven outflow.

The BPT diagram derived for the OC emitting gas shows [N II]/Hα ratios higher than those from NC+OC fluxes, with the highest values associated with the SE and NW nuclear regions. The shift to higher [N II]/Hα in outflowing gas is observed in the SF-driven (e.g. Ho et al. 2014; Newman et al. 2012; Perna et al. 2018a) and AGN-driven (e.g. McElroy et al. 2015; Perna et al. 2017b) outflows. Therefore, it cannot be used to unambiguously distinguish between composite AGN-H II and AGN ionisation. Different mechanisms have been proposed to explain these higher flux ratios: they can be interpreted as due to shock excitation (e.g. Allen et al. 2008), hard ionisation radiation field (assuming that the outflowing material is close to the AGN; e.g. Belfiore et al. 2016), or they could be associated with high-metallicity regions within the outflow (Villar-Martín et al. 2014). The electron temperatures associated with outflowing gas (Table 2) are much lower than those expected for shock excitation (T_e >  5 × 10⁵ K; Osterbrock & Ferland 2006), and may exclude the presence of shock excitation in Mkn 848. On the other hand, metallicity effects cannot be responsible for such high [N II]/Hα ratios (e.g. Curti et al. 2017; Shapley et al. 2015). We therefore suggest that a more plausible explanation for the higher [N II]/Hα in the ejected gas could be given by hard ionisation radiation coming from the central SMBHs.

The figure also shows the [O III]λ5007/Hβ vs. [S II]λ6718,32/Hα flux ratios (central panels) and [O III]λ5007/Hβ vs. [O I]λ6300/Hα flux ratios (bottom panels). These diagrams also allow us to isolate the so-called low-ionisation nuclear emission line regions (LINERs; Heckman 1980), which might be associated with shock excitation (e.g. Allen et al. 2008; but see Belfiore et al. 2016). The NC+OC flux ratios diagnostics confirm AGN-ionised gas in biconical regions in the SE system and H II-like regions in the NW source. The OC diagnostics indicate higher [S II]λ6718,32/Hα and [O I ]λ6300/Hα ratios with respect to the NC+OC values; no clear spatial separation between AGN and shock excitation is found by comparing the different diagrams (see the insets in the panels), however.

We also studied the correlation between the ejected gas kinematics, which are traced by the [O III] W80, and the ionisation state, which is traced by the [N II]/Hα, [S II]/Hα and [O I]/Hα (Dopita & Sutherland 1995; see e.g. Arribas et al. 2014; Mingozzi et al. 2019; Perna et al. 2017a; Rich et al. 2015). We found no clear correlation: all OC flux ratios seem to be independent of the gas kinematics, as expected for photo-ionised gas. Together with the arguments described above, this allows us to prefer AGN hard ionisation radiation over shocks in outflowing gas as a plausible explanation.

To summarise, the BPT diagnostics allow us to determine the nature of the ionised outflows in the AGN-dominated SE system; on the other hand, the NW nucleus has both a strong nuclear SB and an obscured AGN, and it is still not clear if the observed outflow can be unambiguously associated with AGN or SB activity, although the inferred velocities and the higher OC flux ratios seem to suggest that we observe AGN-driven outflows.

9. Outflow energetics

Outflow mass rate, kinetic power, and mass-loading factor μ, which is defined as the ratio between mass outflow rate and SFR, are important ingredients of galaxy evolution models (e.g. Oppenheimer & Davé 2008; Hopkins et al. 2012). These quantities have also been used to distinguish between AGN- and SF-driven outflow (e.g. Brusa et al. 2015b; Harrison et al. 2012). In this section we report the ionised and atomic outflow properties associated with the ejected material observed in the two systems.

Taking advantage of the information so far collected regarding the spatial configuration of the two merging galaxies, the spatial distribution of high-v atomic (Na ID) and ionised ([O III]) gas as well as the absorbing material (traced by Balmer decrement), and the ionisation conditions across the FOV, we show in Fig. 14 a cartoon illustrating the inferred properties of Mkn 848. For the ionised and atomic outflows, our findings are compatible with a simple biconical distribution (e.g. Maiolino et al. 2012).

Fig. 14. Cartoon illustration for Mkn 848. The blue ellipses and black and yellow stars represent the galaxies and the central SMBHs, respectively, with the near edge-on orientation for the SE system and the almost face-on configuration of the NW galaxy. The cyan and grey areas refer to the tidal arms and the dust observed in the HST images. The extended and tenuous yellow area indicates the regions associated with AGN-ionisation: in the SE system, we observe an extended NLR close to the plane of the sky; in the NW, the almost-face configuration suggests that the NLR axis is close to the LOS. The clouds indicate outflows and distinguish between ionised (yellow contour) and atomic (brown contour) gas, as well as between receding (red filled) and approaching (blue filled) material. Clouds visualise outflow directions, phases, and positions inferred from our analysis and do not refer to the (unknown) internal structure of the biconical outflows. The collected information suggest AGN-driven biconical outflows in the two merging galaxies.

Following the arguments presented in Cresci et al. (2015a), we therefore computed the (biconical) outflow mass rates of the ionised components based on the equation

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{H}\beta }^{\mathrm{out}}=8.6\times \frac{L_{41}([\mathrm{H}\beta ])v_{\mathrm{out}}}{N_{\mathrm{e}}{R}_{\mathrm{kpc}}}{M}_{\odot }{\mathrm{yr}}^{-1},\end{array}$$(3)

where L₄₁(Hβ) is the Hβ luminosity associated with the outflow component in units of 10⁴¹ erg s⁻¹ and corrected for extinction, N_e is the electron density, v_out is the outflow velocity, and R_kpc is the radius of the outflowing region in units of kiloparsec (kpc). The coefficient in Eq. (3) was computed assuming T_e  =  2 × 10⁴ K.

For comparison, we also computed the outflow mass rates of the ionised components from the equation presented by Cano-Díaz et al. (2012), who employed the [O III] luminosity:

$$\begin{array}{c}\hfill {\dot{M}}_{[\mathrm{O}\mathrm{III}]}^{\mathrm{out}}=0.48\times \frac{C{L}_{41}([\mathrm{O}\mathrm{III}])v_{\mathrm{out}}}{N_{\mathrm{e}}{R}_{\mathrm{out}}{10}^{\left[\mathrm{O}/\mathrm{H}\right]}}{M}_{\odot }{\mathrm{yr}}^{-1},\end{array}$$(4)

where L₄₁([O III]) is the [O III] luminosity associated with the outflow component in units of 10⁴¹ erg s⁻¹ and corrected for extinction, C is the condensation factor (≈1), N_e is the electron density, and 10^[O/H] is the metallicity in solar units, R_out is the radius of the outflowing region in units of kpc. The coefficient in Eq. (4) was computed considering the dependence of the [O III] emissivity on the electron temperature (2 × 10⁴ K)⁵.

The kinetic and momentum powers are derived from the relations

$$\begin{array}{c}\hfill {\dot{E}}^{\mathrm{kin}}=\frac{1}{2}{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^2=3.2\times {10}^{35}\left(\frac{{\dot{M}}_{\mathrm{out}}}{M_{\odot }{\mathrm{yr}}^{-1}}\right){\left(\frac{v_{\mathrm{out}}}{{\mathrm{km}\mathrm{s}}^{-1}}\right)}^2{\mathrm{erg}\mathrm{s}}^{-1},\end{array}$$(5)

$$\begin{array}{c}\hfill {\dot{P}}^{\mathrm{ion}}={\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}=6.3\times {10}^{30}\left(\frac{{\dot{M}}_{\mathrm{out}}}{M_{\odot }{\mathrm{yr}}^{-1}}\right)\left(\frac{v_{\mathrm{out}}}{{\mathrm{km}\mathrm{s}}^{-1}}\right)\mathrm{dyne}.\end{array}$$(6)

The outflow energetics were derived assuming a simple single radius wind because of the low spatial resolution and the unknown detailed geometrical configuration of the outflow. All the quantities required to derive the outflow properties, as discussed in the previous sections, are reported in Table 3 together with the outflow mass rates, and kinetic and momentum powers. The AGN ionisation does not allow us to derive metallicity measurements for the outflowing gas; therefore we assumed a solar metallicity to derive [O III]-based outflow energetics. We note that these [O III]-based estimates are ≈1 dex lower than those obtained with Hβ, in agreement with previous results (see e.g. Perna et al. 2015a; Carniani et al. 2015). In the following, we therefore refer only to the Hβ-based outflow energetics for the ionised component of the multi-phase outflows.

The contribution to outflow energetics due to neutral ejected material was derived assuming the same geometrical configuration of ionised gas, starting from the equation (Rupke et al. 2002)

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{out}}^{\mathrm{H}}=21\left(\frac{\mathrm{\Omega }}{4\pi }\right)C_f\left(\frac{N_{\mathrm{H}}}{{10}^{21}{\mathrm{cm}}^{-2}}\right)\left(\frac{R_{\mathrm{out}}}{1\mathrm{kpc}}\right)\left(\frac{v_{\mathrm{out}}}{200\mathrm{km}{\mathrm{s}}^{-1}}\right)M_{\odot }{\mathrm{yr}}^{-1},\end{array}$$(7)

where Ω is the solid angle subtended by the wind, C_f is the covering factor, N_H is the column density, R_out is the radius of the neutral outflow, and v_out is the velocity of the ejected Na ID. The related kinetic and momentum powers were derived using Eqs. (5) and (6).

In deriving the atomic outflow energetics for the NW system, we assumed that the neutral outflow is related to the two distinct components travelling at different velocities⁶ (Fig. 11), and that v_out is associated with the maximum velocity V₀₂ of a given component. All other measurements required to compute the outflow energetics and the derived values are reported in Table 3.

Our results suggest that the neutral and ionised outflows in the SE system are likely to have similar masses (Ṁ^out ≈ 2 M_⊙ yr⁻¹) and outflow energetics (${\dot{E}}^{\mathrm{kin}}\approx {10}^{41}$ erg s⁻¹ and Ṗ ≈ 10³³ dyne), while the NW outflow is dominated by the neutral component: the neutral mass outflow rate (Ṁ^out ≈ 10 M_⊙ yr⁻¹), kinetic (${\dot{E}}^{\mathrm{kin}}\approx {10}^{42}$ erg s⁻¹), and momentum (Ṗ ≈ 10³⁵ dyne) powers are one or two dex higher than those associated with ionised gas. Discrepancies between ionised and atomic components have previously been reported for the NW system by Rupke & Veilleux (2013), and have also been found in other nearby systems (Rupke et al. 2017, ∼2 ÷ 3 dex; see also Fiore et al. 2017). We also note that the outflow in the NW system is associated with higher mass rates and powers.

Rupke & Veilleux (2013) detected both ionised and atomic outflows only in the NW system. With respect to their work, we used different techniques to derive outflow mass rates and outflow energetics for the ionised and the atomic components. Rupke & Veilleux (2013) used the Hα emission to trace the ionised outflow, while we preferred to use the Hβ emission, which is not affected by blending problems; they also assumed a very low electron density (10 cm⁻³) to derive the ionised outflow properties, while we measured densities > 1000 cm⁻³ in the NW outflowing gas. Our neutral outflow energetics are instead consistent with those derived in Rupke & Veilleux (2013) within a factor of a few.

In order to investigate the possible origin of the wind, we compared our inferred values of the total outflow kinetic power with the AGN bolometric luminosity and the expected kinetic power ascribed to stellar processes. The latter term was assumed to be proportional to the SFR, and is at most ∼7 × 10⁴¹× (SFR M_⊙ yr⁻¹) erg s⁻¹, following Veilleux et al. (2005). We also derived the mass-loading factors associated with these outflows. SFR measurements were taken from Yuan et al. (2018): SFR^SE = 11 M_⊙ yr⁻¹ and SFR^NW = 76 M_⊙ yr⁻¹. The AGN bolometric luminosities of the two systems are only poorly constrained (Sect. 2); we therefore decided to consider the (extinction-corrected, NC+OC) [O III] luminosity and a standard bolometric correction (3 × 10³; Heckman et al. 2004) to obtain a rough estimate of the bolometric luminosities of the NW and SE AGN: log(L_bol(NW)) ≈ 45.3 erg s⁻¹ and log(L_bol(SE)) ≈ 44.2 erg s⁻¹⁷.

The outflow kinetic power of the NW and SE outflows can be associated with SF-driven winds, assuming a ∼10% coupling between the stellar processes and the observed winds, and a mass-loading factor μ ≈ 0.5. The latter value is consistent with those derived for local luminous IR galaxies (Arribas et al. 2014; Chisholm et al. 2017; Cresci et al. 2017) and high-z star-forming galaxies (Genzel et al. 2014; Newman et al. 2012; Perna et al. 2018a) without evidence of an accreting SMBH.

Similarly, the AGN luminosities in the two systems are high enough to sustain the inferred outflow powers (Table 3): the value of ${\dot{E}}_{\mathrm{kin}}$ is a few 0.1% of the AGN bolometric luminosity for the SE and NW systems, in agreement with models, which predict a reasonable coupling between the energy released by the AGN and the energy required to drive the outflow (e.g. King 2005; Ishibashi et al. 2018; see also Harrison et al. 2018). These outflow energetics therefore do not allow us to unambiguously distinguish between SB- and AGN-driven outflows.

10. Conclusions

We presented a spatially resolved spectroscopic analysis of Mkn 848, which is a complex system of two merging galaxies at z ∼ 0.04. We analysed the publicly available MaNGA data to study the kinematic of warm and cool gas in detail. We distinguished between gravitational motions and outflow signatures and revealed for the first time multi-phase galaxy-wide outflows in the two merging galaxies. In particular, we revealed ionised and neutral ejected gas in the NW and SE systems; we also detected the Na ID emission across the MaNGA FOV, which is also associated with receding outflow gas in the two merging galaxies.

We tested whether the ejected material can be associated with AGN winds or stellar processes, taking advantage of a detailed analysis of the ionised and atomic gas kinematic, ionisation conditions, and outflow energetics. The energetics of the NW and SE outflows did not allow us to unambiguously distinguish between SB- and AGN-driven outflows: the derived mass-loading factors and outflow kinetic powers would require a ∼10% coupling between stellar processes and the winds to explain the outflow measurements; much smaller coupling factors would instead be required in case of AGN winds. Similarly, the BPT diagnostics did not allow us to definitely distinguish between AGN and SB ionisation. Nevertheless, the collected information might suggest an AGN origin of the observed multi-phase outflows: among them are the high-outflow velocities, of the order of 500–1000 km s⁻¹; the BPT flux ratios associated with the outflow components, significantly higher than those of systemic gas; and the absence of any evidence of shock ionisation. Further deeper multi-wavelength observations are required to better constrain the nature of these multi-phase outflows, however.

Our analysis showed that detailed multi-wavelength studies are required to comprehensively characterise the outflows. We found that the neutral outflow component might be related to mass outflow rates and kinetic and momentum powers similar to or even higher than those of the ionised component. This result further emphasises the need for follow-up observations aimed at detecting and characterising the different outflow phases in individual targets.

We finally note that evidence of vigorous and efficient outflows in early- or intermediate-stage merger systems has been collected for several nearby ULIRG/AGN galaxies (e.g. Feruglio et al. 2013, 2015; Rupke & Veilleux 2013, 2015; Saito et al. 2018) and high-z SMG/QSOs (e.g. Perna et al. 2018b and refs therein). These findings suggest that feedback phenomena could be important starting at the early phases of the SB-QSO evolutionary sequence, efficiently depleting the gas reservoirs in the systems (as observed in z >  1 CT QSOs with gas fractions much lower than in normal galaxies with similar properties). This means that these results might be at odds with the scenario in which the important feedback phase follows the SF-dominated phases that are related to the initial merger stages (e.g. Hopkins 2012).

Acknowledgments

We acknowledge the anonymous referee for their constructive comments, which significantly contributed to improving the quality of the paper. GC acknowledges support by INAF/Frontiera through the “Progetti Premiali” funding scheme of the Italian Ministry of Education, University, and Research, and INAF PRIN-SKA 2017 programme 1.05.01.88.04. We acknowledge financial support from INAF under contract PRIN-INAF-2014 (“Windy Black Holes combing Galaxy evolution”). Part of this research was carried out at the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe” during the program “In & Out: What rules the galaxy baryon cycle?” (July 2017). MP thanks M. Dambrosio, who provided a useful and rare insight into 3D projection. 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.

References

All Tables

All Figures

Fig. 1. Chandra 0.5–8 keV cutout (left panel) and HST images with colours in log scale. In the Chandra image, blue contours represent F435W fluxes. In the F435W image, we labelled the two major tidal tails, T1 and T2, that extend over tens of kpc and possibly overlap on the right of the SE nucleus in the figure. A shell-like region (S) is also indicated. In the three-colour cutout, the white boxes show the GMOS FOV from which previous spatially resolved kinematic results have been obtained by Rupke & Veilleux (2013). The white contour represents the S/N = 40 level in the stellar continuum map at ∼5550 Å and is taken as a reference in the next figures. The MaNGA field of view is also shown with white dashed lines.

In the text

Fig. 2. Portion of MANGA integrated spectra of Mkn 848, for the NW and SE regions, with superimposed best-fitting components. Top panels: [O III]+Hβ emission lines from the NW (left panel) and SE (right panel) regions. Central panels: He I emission line and Na ID absorption line system. Bottom panels: [N II], Hα and [S II] lines. The different sets of Gaussian components are shown with solid green curves (NC), dashed and dot-dashed blue lines (OC). Stellar continuum emission has been previously modelled with the pPXF routine and subtracted. The dotted lines mark the wavelengths of Hβ and [O III]λλ4959,5007, in the top panels, He I and Na ID doublet lines in the central panels, and [N II]λ6548, Hα, [N II]λ6581 and [S II] doublet lines in the bottom panels. The systemic redshift is defined on the basis of the best-fit solution for the NC in the NW integrated spectrum. The SE systemic lines show a velocity offset of ∼90 km s−1. The best-fits reproduce all but the Na ID NC in the NW integrated spectrum with high significance.

In the text

Fig. 3. NW (left panel) and SE (right panel) outflow-integrated spectra extracted from the spaxels associated with high-velocity gas in the neutral and/or ionised phases (ΔV02([O III]) <  −550 km s−1 and/or ΔV02(Na ID) <  − 900 km s−1, for the NW system; ΔV98([O III]) >  450 km s−1 and/or ΔV02(Na ID) <  − 750 km s−1, for the SE galaxy; see Sect. 4.2). The shell-like region spaxels are not included in these integrated spectra. The stellar continuum emission has been previously modelled and subtracted from each spaxel. The red curves overplotted on the spectra represent the best-fit solutions that reproduce the line profiles; NC and OC Gaussians are also shown as in Fig. 2. The outflow-integrated spectra extracted from the original MaNGA datacube (hence without pPXF best-fit model subtraction) and the integrated pPXF best-fit spectra are also shown in the figure with solid and dotted grey curves, respectively. For the two grey curves we applied a downward shift for visualisation purposes.

In the text

	Fig. 4. Top-left panel: number of sets of Gaussians we used to model individual spaxel spectra. Top-right and bottom panels: Hβ, [O III], and Hα integrated fluxes obtained from simultaneous best fits (we show only the spaxels with an S/N >  3). Solid contours refer to S/N levels of 10 and 40 in the stellar continuum map (at ∼5550 Å). Black crosses indicate the nuclear regions.
In the text

Fig. 5. High-velocity [O III]λ5007 emission line flux across the FOV obtained by integrating for each spaxel the spectra in the regions 4890 − 5004 Å and 5014 − 5030 Å, which correspond to absolute velocities > 250 km s−1 with respect to the reference frame of the NW nucleus. The map was obtained after subtracting the stellar contribution (using pPXF best-fit spaxel spectra). Solid contours refer to the S/N = 10 and 40 levels in the stellar continuum map. Different spectra, extracted from the regions labelled in the flux maps, are also shown in the insets, with normalised fluxes; blue and red spectra in inset (b) show the most asymmetric line profiles and are extracted from the 3 × 3 pixel regions shown in the central panel with blue and red squares, respectively. Vertical dashed lines indicate the (NW nucleus systemic) wavelength of the [O III]λ5007 line; dotted lines refer to ±250 km s−1 with respect to the NW nucleus systemic velocity.

In the text

Fig. 6. Central panels: Na ID absorption (Abs.) and emission (Em.) line fluxes across the FOV obtained by integrating for each spaxel the spectra in the regions 5860 − 5910 Å and 5890 − 5920 Å, respectively. Negative and positive fluxes were obtained after subtracting the stellar contribution (pPXF best fit); for the absorption fluxes, we also removed the He I emission line (simultaneous best-fit results). Insets: different spectra extracted from the regions labelled in the flux maps, with normalised fluxes: absorption- (emission-) dominated Na ID line profiles are shown in the right (left) panels. Red curves represent the Na ID stellar contribution fitted with pPXF, rescaled for visual purposes. Vertical dashed lines indicate the (NW nucleus systemic) wavelengths of the He I and Na ID doublet lines.

In the text

Fig. 7. Left panel: stellar velocity map derived from pPXF best-fit models. Central panel: [O III] non-parametric peak velocity Vp estimate (an identical map is obtained for the Hα line). Right panel: difference between [O III] peak velocity (Vp) and stellar velocity (Vstar); the most significant discrepancy between the two kinematic tracers is found in regions A and B in the central panel. All velocities are derived taking the narrow-component velocity of the NW -integrated spectrum as a reference.

In the text

Fig. 8. Maps across the FOV of the maximum blue- (left panel) and red-shifted (central panel) velocities for the [O III] best-fit line profile. In particular, the two panels show the differences between maximum velocities (V02 and V98) and stellar velocity (Vstar), as labelled, to reveal outflow signatures. In the V98 − Vstar map the shell region (S) associated with both receding and approaching high-v gas is shown. The right panel shows the [O III] velocity width W80.

In the text

	Fig. 9. Sodium-line velocity maps: the maximum velocity ΔV02 (left panel), highlighting the receding-absorption components, and ΔV98 (central panel), highlighting the approaching-emitting components. The right panel shows the absorption component W80 velocity map.
In the text

	Fig. 10. Balmer decrement f(Hα)/f(Hβ) ratio across the FOV for the total line profiles (left panel) and for the NC (centre) and OC (right panel) components separately. Higher flux ratios correspond to higher dust extinction.
In the text

Fig. 11. Portion of the outflow-integrated spectra of Mkn 848 around the He I and Na ID lines for the NW (left panel) and SE (right panel) regions. The spectra are extracted from the spaxels associated with high-velocity gas in the neutral phase, with ΔV02 <  −900 km s−1, for the NW system, and ΔV02 <  −750 km s−1 for the SE system. The two spectra have a normalised continuum, and the superimposed best-fitting Na ID components in absorption are obtained following Rupke et al. (2005a), assuming a geometrical configuration with partially overlapping sodium atoms (Sect. 6). The He I emission line is instead fitted with a Gaussian profile.

In the text

	Fig. 12. [S II] flux ratio maps derived for the total line profiles (left panel), and distinguishing between systemic (centre) and outflowing gas (right panel). Higher flux ratios correspond to lower electron densities.
In the text

Fig. 13. Top panels: standard diagnostic [O III]/Hβ vs. [N II]/Hα diagrams. The lines drawn in the figures display the curves used to separate H II-like (below the dashed line) from AGN-ionised regions (above the solid line); the locus in between the two curves is associated with composite SF-AGN regions. The solid line is adapted from Kewley et al. 2001; the dashed curve is taken from Kauffman et al. (2003). In the left panel, flux ratios are obtained considering the total line profiles (i.e. NC+OC fluxes); in the right panel, only the OC fluxes have been used. Small colour-coded dots refer to single spaxel measurements associated with the spatial regions shown in the insets (axis scales are the same as in all previous maps). Black spaxels are associated with flux ratios below the curve from Kauffman et al. (2003); purple to red instead show the increasing distance from the same curve. Central panels: standard diagnostic [O III]/Hβ vs. [S II]/Hα diagrams. The lines separates H II-like ionised regions (below the solid line) from AGN-ionised gas (above the dotted and dashed lines) and LINERs (in the right part of the diagrams). The solid and dashed lines are taken from Kewley et al. (2001) and Kewley et al. (2006), respectively. Colour-coded flux ratio dots highlight the separation between AGN, LINERs, and H II regions, and refer to the spatial regions shown in the insets. Bottom panels: standard diagnostic [O III]/Hβ vs. [O I]/Hα diagrams. The colours of the flux ratio dots and separation lines are the same as those in the central panels. In all figures, the largest blue and red dots refer to the flux ratios obtained from the NW and SE outflow-integrated spectra, respectively; only the flux ratios associated with Hβ and [N II] (or [S II], in the central panels; or [O I] in the bottom panels) lines with an S/N >  3 in are shown in the figures. Representative 1σ errors are shown in the bottom corners of the plots.

In the text

Fig. 14. Cartoon illustration for Mkn 848. The blue ellipses and black and yellow stars represent the galaxies and the central SMBHs, respectively, with the near edge-on orientation for the SE system and the almost face-on configuration of the NW galaxy. The cyan and grey areas refer to the tidal arms and the dust observed in the HST images. The extended and tenuous yellow area indicates the regions associated with AGN-ionisation: in the SE system, we observe an extended NLR close to the plane of the sky; in the NW, the almost-face configuration suggests that the NLR axis is close to the LOS. The clouds indicate outflows and distinguish between ionised (yellow contour) and atomic (brown contour) gas, as well as between receding (red filled) and approaching (blue filled) material. Clouds visualise outflow directions, phases, and positions inferred from our analysis and do not refer to the (unknown) internal structure of the biconical outflows. The collected information suggest AGN-driven biconical outflows in the two merging galaxies.

In the text