# The JWST/NIRSpec view of the nuclear region in the prototypical merging galaxy NGC 6240

M. Ceci, G. Cresci, S. Arribas, T. Böker, A. J. Bunker, S. Charlot, K.

Fahrion, I. Lamperti, A. Marconi, M. Perna, et al.

► To cite this version:

M. Ceci, G. Cresci, S. Arribas, T. Böker, A. J. Bunker, et al.. The JWST/NIRSpec view of the nuclear region in the prototypical merging galaxy NGC 6240. Astronomy & Astrophysics - A&A, 2025, 695, 10.1051/0004-6361/202452207. insu-05101691

HAL Id: insu-05101691

https://insu.hal.science/insu-05101691v1

Submitted on 7 Jun 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The JWST/NIRSpec view of the nuclear region in the prototypical merging galaxy NGC 6240

M. Ceci<sup>1,2,*</sup> G. Cresci<sup>2</sup>, S. Arribas<sup>3</sup>, T. Böker<sup>4</sup>, A. J. Bunker<sup>5</sup>, S. Charlot<sup>6</sup>, K. Fahrion<sup>7</sup>, I. Lamperti<sup>1,2</sup>, A. Marconi<sup>1,2</sup>, M. Perna<sup>3</sup>, G. Tozzi<sup>8</sup>, and L. Ulivi<sup>1,2,9</sup>

Received 11 September 2024 / Accepted 3 January 2025

## ABSTRACT

Merger events are thought to be an important phase in the assembly of massive galaxies. At the same time, active galactic nuclei (AGN) play a fundamental role in the evolution of their star formation histories. Both phenomena can be observed at work in NGC 6240, a local prototypical merger classified as an ultraluminous infrared galaxy (ULIRG) thanks to its elevated infrared luminosity. Interestingly, NGC 6240 hosts two AGN separated by 1.5" (~735 pc), detected in both X-ray and radio band. Taking advantage of the unprecedented sensitivity and wavelength coverage provided by the integral field unit (IFU) of the NIRSpec instrument on board JWST, we observed the nuclear region of NGC 6240 in a field of view of 3.7" × 3.7" (1.9 × 1.9 kpc<fnref n="2" />) in order to investigate gas kinematics and interstellar medium (ISM) properties with a high spatial resolution of ~0.1" (or ~50 pc). We characterized the 2D stellar kinematics, separated the different gas kinematic components through multi-Gaussian fitting, and studied the excitation properties of the ISM from the near-infrared diagnostic diagram based on the H₂ 1-0 S(1)/Bry and [Fe II] λ1.257 µm/Paβ lines ratios. We isolated the ionization cones of the two nuclei and detected coronal line emission from both of them. Using H₂ line ratios, we found that the molecular hydrogen gas is excited mostly by thermal processes. We computed a hot molecular gas mass of 1.3 × 10⁵ M☉ and an ionized gas mass in the range of 10⁵–10⁷ M☉, depending on assumptions. We studied with unprecedented spatial resolution and sensitivity the kinematics of the molecular and ionized gas phases, and we revealed the complex structure of the molecular gas and found a blueshifted outflow near the southern nucleus, together with filaments connecting a highly redshifted H₂ cloud with the two nuclei. We speculate on the possible nature of this H₂ cloud and propose two possible scenarios: outflowing gas or a tidal cloud falling onto the nuclei.

**Key words.** galaxies: active – galaxies: interactions – galaxies: nuclei – quasars: emission lines – quasars: supermassive black holes – infrared: galaxies

## 1. Introduction

Ultraluminous infrared galaxies (ULIRG) are among the most luminous objects in the local Universe, and they represent a fundamental phase in the evolution of active galactic nuclei (AGN). This extremely luminous episode, due to a combination of a starburst and/or AGN activity in a dust-enshrouded environment, occurs probably during a major merger of gas-rich galaxies (Sanders et al. 1988; Sanders & Mirabel 1996; Lípari et al. 2003; Hopkins et al. 2006), when gas and dust are both channeled into the nuclear regions of the progenitors (U et al. 2012; Lonsdale et al. 2006). In this scenario, ULIRGs are seen as a rapid growth phase when the merger initiates an intense period of star formation, fueling the growth of a supermassive black hole (SMBH) and the eventual AGN activity. Throughout the merger process, powerful outflows generated by both starburst and AGN activity shape the galactic surrounding, influencing subsequent

star formation and black hole accretion through feedback mechanisms (e.g., Silk 2013, Alexander & Hickox 2012).

In this paper, we focus on NGC 6240. Despite having an IR luminosity of 6.3 × 10¹<fnref n="1" /> L☉, which is just below the canonical ULIRG luminosity boundary of 10¹<fnref n="2" /> L☉, it is frequently associated with the ULIRG category and stands out as one of the closest and most extensively examined objects within this category (e.g., Genzel et al. 1998). At z ~ 0.02, NGC 6240 is considered the prototypical merger galaxy in the local Universe. This means it is observed during the short-lived and therefore rarely observed transient phase between the first encounter and the final coalescence of two gas-rich spiral galaxies (Tacconi et al. 1999; Engel et al. 2010). It hosts two X-ray confirmed AGN that are positioned approximately 735 pc apart (Komossa et al. 2003; Max et al. 2007) and driving kiloparsec-scale superwinds detected across various wavelengths: radio (Baan et al. 2007), far-IR (Veilleux et al. 2013), near-infrared (NIR; van der Werf et al. 1993, Max et al. 2005), optical (Heckman et al. 1990; Veilleux et al. 2003), and X-ray

<footnote n="1">¹ Università di Firenze, Dipartimento di Fisica e Astronomia, via G. Sansone 1, 50019 Sesto F.no, Firenze, Italy</footnote>

<footnote n="2">² INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="3">³ Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir Km. 4, 28850 Torrejón de Ardoz, Madrid, Spain</footnote>

<footnote n="4">⁴ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="5">⁵ Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="6">⁶ Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="7">⁷ Department of Astrophysics, University of Vienna, Türkenschanzstraße 17, 1180 Wien, Austria</footnote>

<footnote n="8">⁸ Max-Planck-Institut für Extraterrestrische Physik (MPE), Giessenbachstraße 1, D-85748 Garching, Germany</footnote>

<footnote n="9">⁹ University of Trento, Via Sommarive 14, I-38123 Trento, Italy</footnote>

<footnote>* Corresponding author; matteo.ceci@unifi.it</footnote>
(Nardini et al. 2013; Wang et al. 2014). Despite both being active, the two nuclei might not exert enough influence or may be too obscured to ionize the interstellar medium (ISM) significantly. Indeed, star formation likely contributes 76–80% of the IR luminosity (Armus et al. 2006; Nardini et al. 2009).

NGC 6240 exhibits one of the most intense NIR $H_2$ emission observed in any galaxy, with a luminosity ($L_{H_2}$) of approximately $10^8 L_\odot$ (Joseph et al. 1984). Cicone et al. (2018) found a molecular outflow extending more than 10 kpc along the E–W direction and peaking between the two AGN rather than on either of the two. Numerous studies have charted the $H_2$ line emission within the inner (~1 kpc) region (e.g., Tecza et al. 2000; Max et al. 2005; Ilha et al. 2016), indicating that shock excitation is the predominant mechanism. Treister et al. (2020) and Fyhrie et al. (2021) argue that their observations are consistent with a tidal bridge between the two nuclei and an observed high velocity redshifted gas. Gas that is subject to gravitational torques, such as that in a tidal bridge in the nuclear region, will likely fall into the nuclear regions by the time of final coalescence (Souchay et al. 2013). However, it is also possible that the structure will dissipate before this streaming occurs, given the high velocity dispersion relative to the vertical extent of the gas.

One of the most peculiar features of NGC 6240 concerns the positions of the two galaxy nuclei. The distance between them appears consistent across different wavelengths: the two AGN are separated by approximately 1.5" (~700 pc) in X-ray (Komossa et al. 2003), radio (Gallimore & Beswick 2004), and at 3–5 µm (Risaliti et al. 2006). The first high-resolution view of the NIR nuclei by Max et al. (2005) revealed the impact of dust on the shapes and positions of the nuclei from optical to 2.2 µm, with Max et al. (2007), demonstrating the consistency between the radio/X-ray positions and those at high resolution in the NIR (separation ~1.5"). However, the bright spots seen in the NIR and optical, separated by about 1.8" (~900 pc; e.g., Gerssen et al. 2004), are significantly affected by dust extinction. As a result, these spots do not correspond to the true positions of the nuclei, which only become visible at wavelengths longer than 2–3 µm, where the dust obscuration is less pronounced. Kollatschny et al. (2020) propose the presence of a third SMBH close to the southern one based on narrow field mode (NFM) MUSE data. Fabbiano et al. (2020) investigated the existence of this putative third nucleus through a reanalysis of the cumulative ACIS Chandra data but did not find any strong evidence for it, so the question is still under discussion.

This paper presents new observations of the nuclear region of NGC 6240 obtained with the integral field unit (IFU) of the Near Infrared Spectrograph (NIRSpec) instrument on board the James Webb Space Telescope (JWST; Jakobsen et al. 2022; Böker et al. 2022) and collected as part of the JWST/NIRSpec IFS Guaranteed Time Observations (GTO) survey Resolved structure and kinematics of the nuclear regions of nearby galaxies (program lead: Torsten Böker). Other galaxies studied in this program include NGC 4654 (Fahrion et al. 2024), Arp 220 (Perna et al. 2024; Ulivi et al. 2025), and Mrk 231 (Ulivi et al. in prep.). In collaboration with the Mid-Infrared Instrument (MIRI) GTO team, this program offers a complete NIR and mid-IR perspective on the nuclear region of NGC 6240 (see the MIRI NGC 6240 paper by Hermosa Muñoz et al. 2025).

This manuscript is organized as follows. Section 2 presents the NIRSpec observations and data reduction. Section 3 presents the spectral analysis using a multi-Gaussian fitting procedure. In Section 4, we discuss our results through the characterization of the ISM properties and the spatially resolved kinematics of the ionized and molecular gas in the central region of NGC 6240.

Finally, Section 5 summarizes our conclusions. Throughout this work, we assume $\Omega_m = 0.286$ and $H_0 = 69.9 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Bennett et al. 2014).

## 2. Observation and data reduction

### 2.1. JWST/NIRSpec data

NGC 6240 was observed with the JWST/NIRSpec instrument in IFU mode as part of the GTO programme 1265 (PI: A. Alonso-Herrero), on 17th July 2023. Centered on an intermediate position between the two nuclei of NGC 6240 (coordinates from far-IR continuum emission; see Scoville et al. 2015), the dataset utilizes the following high-resolution ($R \sim 2700$) configurations, with their respective spectral ranges: [0.97 µm, 1.89 µm] for G140H/F100LP (cube1 hereafter), [1.66 µm, 3.17 µm] for G235H/F170LP (cube2) and [2.87 µm, 5.27 µm] for G395H/F290LP (cube3). The observations were taken with an NRSIRS2RAPID readout pattern with 15 groups per integration, and one integration per exposure, using a four-point medium-cycling dithering. Every configuration has a total integration time of 1867 s. However, this work focuses solely on the first two configurations, as our aim is to study the kinematics of the emitting gas, with its strongest emission lines below 2.4 µm. In this wavelength range, D'Eugenio et al. (2024) estimated a FWHM Point Spread Function (PSF) of ~0.1".

We reduced the data following the same procedure as described in Perna et al. (2023); here we summarize the most relevant details. We used the STScI pipeline v1.8.2 with CRDS context 1063. A patch was incorporated to address specific pipeline problems that affect this specific version of the pipeline; 1/f noise and outliers were corrected following the approaches described in Perna et al. (2023) and D'Eugenio et al. (2024), respectively. The combination of exposures for the four dither positions utilized a drizzle weighting technique, as detailed in Law et al. 2023, which enabled sub-sampling of detector pixels. This process yielded cube spaxels measuring 0.05" each, approximately corresponding to 25 pc/spaxel and covering a field of view (FoV) of 3.7" × 3.7" (1.9 × 1.9 kpc²). The reduced datacubes showed the so-called "wiggles", sinusoidal modulations due to the PSF undersampling by the NIRSpec detectors, affecting mostly the brightest spaxels in the nuclear regions. We corrected them using the same procedure as described in Perna et al. (2023), applying it on the two nuclear regions separately.

For NGC 6240 we measured a redshift of $z = 0.02474 \pm 0.00011(D = 106 \text{ Mpc})$ from a spectrum extracted from an aperture of 30-spaxel radius, centered on the NIRSpec FoV. It was computed as the centroid of the brightest molecular line, that is $H_2$ 1–0 S(1), modeled using 2 Gaussian profiles. This value is compatible with $z = 0.0245$ by Müller-Sánchez et al. (2018). We will use our measurement in the rest of this paper.

To correct the datacubes for astrometry, we assume that the peak of $P\alpha$ in the southern nucleus to be coincident with the southern radio source observed by Gallimore & Beswick (2004), because we expect the emission line to peak at the position of the AGN. We applied astrometric corrections to the three NIRSpec cubes, resulting in shifts of $\Delta RA = -0.2252''$ and $\Delta Dec = +0.1750''$. These offsets are consistent with the expected accuracy of JWST pointing in the absence of a dedicated target acquisition procedure (Böker et al. 2023).

To illustrate the quality of the spectral information contained in the NIRSpec datacubes, we show respectively in Figs. 1 and 2 spectra of the northern and southern nucleus integrated
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Spectrum of the northern nucleus, integrated over an aperture of 5-spaxel radius ($r = 0.25''$ corresponding to ~130 pc). The blue curve is the spectrum from G140H/F100LP (cube1), the orange one from G235H/F170LP (cube2), and the red one from G395H/F290LP (cube3). Each cube spectrum shows a discontinuity in the middle of the band, which is due to the gap between the NIRSpec detectors. We indicated the brightest emission lines with vertical dashed lines (in blue those of molecular hydrogen, in black the others) and the more noticeable features, as the CO absorption bands, polycyclic aromatic hydrocarbons (PAH) and aliphatic emission, and CO$_2$ ice absorption.</figcaption></figure>

over an aperture of radius 5 spaxels ($r = 0.25''$ corresponding to ~130 pc), with labels of the main emission lines and features. The discontinuities in the middle of the three bands are due to the gap between the NIRSpec detectors. The overlapping wavelength ranges at ~1.85 µm and ~3 µm show a good match between the different spectra, demonstrating a consistent flux calibration.

## 2.2. Narrow-Field Mode MUSE data

In order to compute electron density from the optical [S II] doublet (see Section 4.2.5), we also used data obtained and reduced by Kollatschny et al. (2020)<fnref n="1" /> on April 22, 2018, as part of the commissioning run of the MUSE (Multi Unit Spectroscopic Explorer) instrument in its Narrow-Field Mode (NFM) using the four-laser adaptive optics system of ESO's Very Large Telescope (VLT) unit telescope four ("Yes pun"). They took four 500 s science exposures centered on NGC 6240 and another 500 s exposure on an offset sky field. The NFM covers a field of 7.5" × 7.5" on the sky. Between each on-target exposure, they rotated the FoV by 90 degrees. Spatial offsets of about 1.5" provided significantly larger coverage of ~11.5" × 11.5" with some gaps at the field edges. The wavelength range covered is 4725–9350 Å, with a spectral resolution of approximately 2.5 Å. The spectra are sampled at 1.25 Å in the dispersion direction and 0.0253" in the spatial direction. However, there is a gap in the spectrum between 5780 and 6050 Å due to the notch filter to

exclude the laser-guide sodium emission. The data reduction is described in detail in Kollatschny et al. (2020). Similar to NIR-Spec data, we used peaks of $H\alpha$ to correct MUSE observations for astrometry. Corrections result in shifts of $\Delta RA = +1.3944''$ and $\Delta Dec = -1.7892''$.

## 3. Analysis

In this section we describe the methods used to analyze the kinematics and properties of the emitting gas, following Marasco et al. (2020), Tozzi et al. (2021), and Cresci et al. (2023). We divided the fitting procedure into three steps: firstly, we subtracted the stellar contribution from the cubes; then we applied a multiple-component fitting procedure, dividing the gas species into two kinematically independent families; at the end we separated the emission of the narrow component (for example, disk and/or merger debris) from the broad component (outflow and/or turbulent gas). These three steps are described in detail in the following sections.

We applied two consecutive spectral fittings: firstly, we fit the full datacubes aiming to isolate and remove stellar and AGN continuum emission (see Section 3.1); secondly, we modeled the gas emission lines in the continuum-subtracted cubes (see Section 3.2) via a refined modeling with more freedom on the choice of line profiles and model constraints.

### 3.1. Modeling of the stellar continuum

We applied the penalized PiXel-Fitting (pPXF; see Cappellari & Emsellem 2004) package to model the continuum and gas emission, convolving a linear combination of stellar templates

<footnote n="1">¹ The reduced datacubes are available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/633/A79</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Spectrum of the southern nucleus integrated over an aperture of 5-spaxel radius ($r = 0.25''$ corresponding to ~130 pc). See Fig. 1 for further details.</figcaption></figure>

with a Gaussian velocity distribution. First, to increase the signal to noise ratio (S/N) of our datacubes, we used the Voronoi binning algorithm by Cappellari & Copin (2003). The applied S/N threshold per wavelength channel for each bin was 5 for both cube1 and cube2. Regarding the stellar contribution, we adopted the single stellar population models of Maraston & Strömbäck (2011), which have a resolution of $R = 20\,000$ and a spectral range of [0.1 µm, 2.5 µm]. To take into account also the continuum from AGN, we added to the model a 3-degree polynomial. We also fit all the emission lines listed in Table 1. We built a model for the emission lines, with a variable number of Gaussian components (from 1 to 3, because along the line of sight, we can have different kinematic features) and we added it to the previous stellar model, thus obtaining our total model. The kinematic parameters ($v$ and $\sigma$) are the same for all emission lines. We fit each spectrum three times (with 1, 2, or 3 Gaussian components), and then applied a Kolmogorov-Smirnov (KS) test on the residuals to choose the optimal number of Gaussian components in each spaxel, with a threshold $p$-value of 0.5 (Marasco et al. 2020). In practice, we test whether adding an extra component ($n+1$) results in residuals that are statistically different from those of a model with $n$ components, starting with $n=1$. If there is a significant difference, we increase $n$ by 1 and repeat the test. If not, we consider $n$ to be the optimal number of components for that specific spaxel. This method is independent of any particular likelihood estimator and focuses solely on the relative improvements that an additional component brings to the model. The tolerance for accepting a new component can be adjusted by changing the $p$-value threshold of the KS test.

Examples of fit stellar component are shown in Fig. 3, where we plot the fit cube2 spectra of a spaxel centered on the southern nucleus (upper panel) and another located in the region between the two nuclei (lower panel). In both, we show data (in black), the total model of emission lines and stellar component (in red), the model of only stars (in green), and the residuals between data and total model (in blue). The vertical lines mark all features

modeled with pPXF procedure, while the gray bands show the spectral regions we masked during the fitting procedure, corresponding to the gap between the detectors and covering a few emission lines not relevant to this work. We focused on fitting only the most relevant emission lines to ensure a reasonable computational time. We note that the detector gap causes many wavelength channels to be lost and thus prevents a complete fit of the CO band-head absorption features longward of 2.4 µm. In addition, an imperfect wiggles correction has probably spoiled these features. To be more conservative, we decided to include in the fit only the first band-head (~2.35 µm in the observed frame), that is the one furthest from the spectral discontinuity due to the gap between the NIRSpec detectors, the one with fewer residuals by wiggles correction (see Fig. 3). The main differences between the spectra of these two regions are the following. In the nucleus, the continuum flux is about one order of magnitude higher, and the predominant emission line is $Pa\alpha$, while in the region between the nuclei the molecular hydrogen emission lines are the strongest. The somewhat higher residuals in the nuclear spaxel are likely due to an imperfect wiggles correction. The cube1 spectra of the same spaxels, which have lower S/N, are shown in Appendix A.

We plot in Fig. 4 the velocity and velocity dispersion maps of the stellar component, derived from cube2 with pPXF fitting. For these and the other maps in the rest of the paper, north is up and east is left. The velocities are calculated with respect to the systemic one, derived from our fiducial redshift (see Section 2). From the stellar velocity map we can distinguish velocity gradients likely associated with the rotating stellar disks of the two nuclei. The major axis of the southern disk has a position angle$^2$ (PA) of 321.4° ± 6.4°, estimated by fitting the star velocity map with the FIT KINEMATIC PA code (Cappellari et al. 2007, Krajnović et al. 2011). The disk has positive velocities –

<footnote>$^2$ The position angle is measured counterclockwise from the north direction on the approaching part of the disk.</footnote>
<table><caption>**Table 1.** List of all emission lines used during the fitting procedure in NIRSpec datacubes.</caption><thead><tr><th>Line</th><th>λ<sub>vac</sub> (Å)</th><th>Group</th><th>Line</th><th>λ<sub>vac</sub> (Å)</th><th>Group</th></tr></thead><tbody><tr><td>[C I]</td><td>9827</td><td>1</td><td>H<sub>2</sub> 1–0 S(6)</td><td>17880</td><td></td></tr><tr><td>[C I]</td><td>9850</td><td>1</td><td>H<sub>2</sub> 6–4 O(4)</td><td>17960</td><td></td></tr><tr><td>Paδ</td><td>10052</td><td></td><td>[Fe II]</td><td>18099</td><td></td></tr><tr><td>He I</td><td>10314</td><td></td><td>H<sub>2</sub> 1–0 S(5)</td><td>18358</td><td></td></tr><tr><td>[Fe III]</td><td>10806</td><td></td><td>Paα</td><td>18756</td><td>1</td></tr><tr><td>He I</td><td>10832</td><td></td><td>H<sub>2</sub> 1–0 S(4)</td><td>18920</td><td>2</td></tr><tr><td>H<sub>2</sub> 2–0 S(5)</td><td>10851</td><td></td><td>Brδ</td><td>19451</td><td>1</td></tr><tr><td>Paγ</td><td>10941</td><td>1</td><td>H<sub>2</sub> 1–0 S(3)</td><td>19576</td><td>2</td></tr><tr><td>[P II]</td><td>11886</td><td>1</td><td>[Si VI]</td><td>19625</td><td>1</td></tr><tr><td>[Fe II]</td><td>12567</td><td>1</td><td>H<sub>2</sub> 1–0 S(2)</td><td>20338</td><td>2</td></tr><tr><td>Paβ</td><td>12821</td><td>1</td><td>H<sub>2</sub> 2–1 S(3)</td><td>20735</td><td></td></tr><tr><td>[Fe II]</td><td>15335</td><td></td><td>H<sub>2</sub> 1–0 S(1)</td><td>21218</td><td>2</td></tr><tr><td>[Fe II]</td><td>16440</td><td>1</td><td>Bry</td><td>21661</td><td>1</td></tr><tr><td>[Fe II]</td><td>16773</td><td>1</td><td>H<sub>2</sub> 1–0 S(0)</td><td>22235</td><td>2</td></tr><tr><td>H<sub>2</sub> 1–0 S(9)</td><td>16877</td><td>2</td><td>H<sub>2</sub> 2–1 S(1)</td><td>22477</td><td>2</td></tr><tr><td>H<sub>2</sub> 1–0 S(7)</td><td>17480</td><td>2</td><td>H<sub>2</sub> 1–0 Q(5)</td><td>24548</td><td></td></tr></tbody></table>

**Notes.** See Section 3.1 for details of the pPXF modeling. The “group” label defines the two line species used in the fitting procedure part (see Section 3.2).

which means stars moving away from the observer – reaching a measured value of +350 km s⁻¹ in the south-east side, while the side approaching the observer is the north-west one, with a maximum of −250 km s⁻¹ with respect to the systemic velocity. The northern nucleus has more moderate velocities: with PA of 218.6° ± 6.4°, it has the redder part in the NE direction and the bluer one in the SW direction. Respectively, the velocities peaks are +250 km s⁻¹ and −150 km s⁻¹. In addition, the southern disk shows asymmetric velocity curve along the major axis, and in both disks the minor axes seem to be curved. These departures from regular rotation are likely related to gravitational interactions associated to the merging process.

Also the velocity dispersion map is not simply compatible with the presence of two disks, for which we would expect to observe the highest dispersion in the rotation center. Indeed, we find the lowest velocity dispersion of 100–150 km s⁻¹ south of the southern nucleus and east of the northern nucleus, while the highest values of about 350–400 km s⁻¹ are found between the nuclei and in a region to the west of the southern nucleus. These anomalously high dispersions is expected in a merger system. Indeed, these features could be attributed to emission from stars recently formed due to the interaction, which is superimposed on the light from the nucleus and leads to an overestimate of the velocity dispersion. These maps are derived from the cube2 data, but similar and compatible results are found using cube1.

The stellar kinematics of the inner region of NGC 6240 have already been investigated with NIR IFS data from VLT/SINFONI and Gemini (e.g., Engel et al. 2010; Ilha et al. 2016). Fitting the CO absorptions band-heads at ~2.3 µm, Engel et al. (2010) and Ilha et al. (2016) found a velocity amplitude of ~600 km s⁻¹ for the southern nucleus and ~350 km s⁻¹ for the northern one, comparable to our findings of 600 km/s and 400 km/s for the northern and southern nuclei, respectively. Their velocity dispersion maps are similar to ours with values in a range of ~150 km s⁻¹ and ~400 km s⁻¹ in both papers. Regarding the PA of the stellar disks, Engel et al. (2010) derived 329° ± 2° and 229° ± 2° for the southern and northern nucleus,

respectively. Instead, Ilha et al. (2016) estimated PA ~ 330° and PA ~ 225°. Both estimations are comparable with ours within the 2-sigma uncertainties.

Having the stellar model for every spaxel in hand, we subtracted it from the original spaxel spectra, in order to obtain cubes containing only emission lines from the ISM of the galaxy. More precisely, we subtracted the stellar emission from the original unbinned datacube by rescaling the fit stellar model (constant within each Voronoi bin) to the original unbinned spectrum of each spaxel before subtracting it.

## 3.2. Multi-Gaussian emission line fit

After the subtraction of the stellar component, we have datacubes containing just the emission of the ionized, neutral, and molecular gas. In Section 3.1, we used Voronoi binning to combine spaxels and trade spatial resolution to increase S/N. For our emission line analysis, we instead choose to smooth each spatial plane of both our cubes with a Gaussian kernel of size σ<sub>smooth</sub> = 1 pixel, corresponding to 0.05". This kernel is smaller than the PSF (on the order of 0.1", D'Eugenio et al. 2024), and degrades the resulting spatial resolution only minimally.

At this point, we performed the emission line fitting spaxel-by-spaxel. The spectra in the continuum-subtracted datacube were refit for the following reasons: although we previously performed the emission line fitting using the pPXF package in Section 3.1, we are now applying a more refined fitting process with additional constraints and line profile options. We adopted again the strategy of a multiple components fitting procedure, to better reproduce the more asymmetric line profiles. In addition, we separated emission lines into two families, in order to probe potential kinematic differences between the molecular gas (traced by the roto-vibrational transitions of H₂), and the ionized and neutral gas. The advantage of this method is that it overcomes the low S/N of some individual lines by fitting all lines within the same group together. The main disadvantages are that having a relatively high fraction of lines with low S/N may increase the scatter in the results and that we are fitting a large number of lines together, which will be computationally heavy with more components. The obvious work-around would be to fit the two groups in separate runs, but, as we can see in Fig. 3, some spectra contain blended lines with components from different line species that have to be fit simultaneously.

We performed the fit for cube1 including the molecular family of H₂ 1–0 S(7) and H₂ 1–0 S(9), and the other family of [C I]λ0.983 µm, [C I]λ0.985 µm, [P II]λ1.189 µm, [Fe II]λ1.257 µm, [Fe II]λ1.644 µm, Paβ and Paγ. For cube2, we instead fit the group of molecular lines H₂ 1–0 S(2), H₂ 1–0 S(3), H₂ 1–0 S(1), H₂ 1–0 S(4), H₂ 1–0 S(0) and H₂ 2–1 S(1), and the one of [Si VI], Brδ, Paα and Bry. These emission lines are listed in Table 1, where the separation in the two kinematically independent groups of both datacubes is indicated. For all emission lines in either group, the fit parameters for the Gaussian profile are the flux, the velocity v with respect to the systemic velocity, and the velocity dispersion σ. Within each line family and for every spaxel, the fit velocity and velocity dispersion were forced to be identical, while only the line flux was allowed to vary freely between lines. We added also an additive polynomial for cube1 and cube2, respectively of degree 2 and 3, because, there was some residual continuum contribution in the nuclear spatial regions affected by wiggles. We carried out the fit using from one to four Gaussian components, to reproduce all the kinematic components along the line of sight of every spaxel,
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Examples of continuum and emission lines fitting in a spaxel from the southern nucleus (top panel) and one from the region between the nuclei (bottom panel) with three Gaussian components for each line. The data, the emission lines model, the stellar model, and the residuals are shown in black, red, green, and blue respectively. Gray bands indicate spectral regions excluded from the fit. Dashed vertical lines mark the emission lines modeled with pPXF fitting (listed in Table 1).</figcaption></figure>

using for each parameter physical and reasonable boundary conditions<fnref n="3" />.

After the fitting procedure, we have to select the optimal number of components in every spaxel (as done in Section 3.1). We performed again the KS test, changing this time the *p*-value. We increased it to 0.9 to reproduce the line profiles with a higher accuracy. As expected from the chaotic kinematics, 3–4 Gaussian components were necessary in the majority of the spaxels to reproduce the highly asymmetric line profiles. Moreover, there is not any particular shape in the spatial distribution of individual Gaussian components that can be

used to point out some kinematic structure, such as an outflow cone.

After modeling the emission line gas, we aim to decompose the line profiles. This is done by separating the narrow and broad components based on some reasonable criteria.

### 3.3. Decomposition of narrow and broad components

In this section we describe the method used to distinguish between different kinematic components. Achieving a detailed and precise identification of each component's contribution to the overall line emission spaxel-by-spaxel proves challenging, given the complexity of the velocity structure in the FoV. Hence,

<footnote n="3">³ $|v| \le 1000 \text{ km s}^{-1}$ and $\sigma$ greater or equal than the NIRSpec spectral resolution.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Velocity (left panel) and velocity dispersion (right panel) maps of the stellar component derived from cube2 with the pXRF fitting. White contours tracing continuum are overplotted as a spatial reference (for the southern nucleus we used 19% of the 1.2 µm continuum maximum and for the northern nucleus 6% of the 2.3 µm continuum maximum). Crosses and the ellipse indicate respectively the nuclei and the high-velocity molecular hydrogen blob positions (see Section 4.3.1). From the velocity map, we can distinguish the two stellar rotating disks in the proximity of the two nuclei.</figcaption></figure>

we executed a decomposition between broader components possibly linked to outflowing gas, which we refer to as “broad” components, and the others we simply call “narrow”. The criterion we applied is as follows: we defined the wings of the total line profile as the tails exceeding a predefined velocity threshold with respect to the velocity peak of the total line profile. Then, we calculated for each Gaussian component the fraction of flux within these wings relative to the total flux of the Gaussian component. If this fraction exceeds a fixed threshold, the component is classified as broad. Otherwise, it is classified as narrow. The parameters of this method are the velocity threshold that defines the wings of the line, and the fraction of the total flux of the component. We used $v_{\text{wings}} = 200 \text{ km s}^{-1}$ and $f_{\text{wings}} = 0.5$ for both datacubes. This means that a component is classified as ‘broad’ if more than 50% of its flux is at velocities >200 km s$^{-1}$ w.r.t. the peak of the total line profile. Our classification criterion relies on the fact that outflows typically generate broader and higher-velocity line emission (Marasco et al. 2020, Tozzi et al. 2021), contrasting with the other potential components, which usually produce narrow line emission near the galaxy’s systemic velocity.

We present some examples of decomposition in Fig. 5. In every row we show the spectrum of the spaxel marked with a light green point in the first panel on the left, that is the image of the Paα flux in logarithmic scale. In the other panels, we plot the continuum-subtracted spectrum, the total model, the model of the narrow components, and the model of the broad components (in black, blue, green, and orange, respectively) of the brightest emission lines. The narrow and broad sub-components are grouped to better analyze the global properties of the low- and high-velocity gas. The vertical lines represent line wavelength corresponding to the galaxy systemic velocity. As we can see, the shape and relative strength of the two components vary substantially across the FoV. We also find that the ionized gas always has a broad component, while the molecular gas sometimes has only narrow components. As can be seen in Fig. 5 the Paα is sometimes blended with He I. We decided not to fit this latter because it is negligible and when 3–4 Gaussian components are

used, some He I components are classified as Paα outflow (and vice versa).

### 3.4. Nuclei positions

As shown by Max et al. (2005), we also observed that the apparent location of the nuclei shifts as longer wavelengths enable us to peer deeper into the dust to reveal the true continuum centers. Using adaptive optics facilities from the ground, they concluded that in the K' band the AGN in the southern nucleus is still enshrouded in dust and they tentatively identified the northern nucleus position. Our longer wavelength JWST/NIRSpec data enable us to identify the true position of both AGN.

In each column of Fig. 6, we show the continuum images of the datacubes, zooming in on the nuclei positions. White contours of the continuum are overlaid as spatial references, using the 1.2 µm continuum for the southern nucleus and the 2.3 µm continuum for the northern nucleus. The northern nucleus is shown in the top panels and the southern nucleus in the bottom. From left to right, the columns present the continuum at 1.2 µm, 2.3 µm, and 5.0 µm (rest frame) obtained with a window of ±0.01 µm. We can easily observe that the continuum emission peaks shift across the NIR spectral range. Looking at the northern nucleus, we show that the peak position moves from the northeastern position – observable in the optical band (Gerssen et al. 2004) – to the southwestern one, which matches the radio position (Gallimore & Beswick 2004). Regarding the southern nucleus, in the visible band two sub-components are present in the nuclear region. As we move to redder wavelengths, the southern component fades, while the northern component shifts slightly to the north until it reaches the position of the radio nucleus at 5.0 µm. We thus confirm the variation of the position of the continuum peaks through infrared wavelengths, as also reported in previous works (e.g., Gerssen et al. 2004; Max et al. 2005, 2007; see also Scoville et al. 1998 and Perna et al. 2024 for a similar behavior in Arp 220). We verified that the continuum peaks at the longest NIR wavelengths correspond to the positions of X-ray and radio sources, confirming the results of
<figure><img src="image_6.png" /><figcaption>Fig. 5. Examples of the decomposition of narrow and broad components at the position of some spaxels. The spaxels are shown with a green point in the left panels (which are Paα maps in logarithmic scale). Black crosses and the black ellipse indicate respectively the nuclei and the high-velocity molecular hydrogen blob positions. The other panels show the results of the decomposition for the Paβ, [Fe II]λ1.644 µm, Paα. H₂ 1–0 S(4), and H₂ 1–0 S(1). Black curve shows the observed emission lines after the subtraction of stellar continuum, the blue curve shows the total result of the fit, the green and orange curves mark the narrow and broad components, respectively. The systemic velocity is represented by a gray dashed line. The wavelengths on the x-axis are in the rest-frame. The y-axis units are erg s⁻¹ cm⁻² Å⁻¹.</figcaption></figure>

Max et al. (2007) and the validity of our initial choice to use the Paα peak to mark the AGN positions.

### 3.5. Extinction correction

Following Ramos Almeida et al. (2009), we estimated the extinction using the Paβ/Paα and the [Fe II]λ1.257 µm/[Fe II]λ1.644 µm line ratios. In Fig. 7, we plot the observed Paschen and [Fe II] line ratios and then we overplot inferred values for fixed $A_V$, assuming a range of physically reasonable values of gas density and temperature. We

used the Cardelli et al. (1989) extinction law, valid in the range $0.9 \mu m \le \lambda \le 3.3 \mu m$,

$$A(\lambda) = A_V (a(\lambda) + b(\lambda)/R_V),$$

where

$$a(\lambda) = 0.574 \lambda^{-1.61},$$

$$b(\lambda) = 0.527 \lambda^{-1.61},$$

and $R_V = 3.1$, which is the standard value for the diffuse ISM. As shown in the plot, while iron ratios are not very
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Zoomed-in continuum maps of the northern (top) and southern (bottom) nuclear region for different rest-frame wavelengths, as indicated by the axis. Crosses indicate the positions of the nuclei. White contours of the continuum flux are overplotted as a spatial reference (for the southern nucleus we used 19% of the 1.2 µm continuum maximum and for the northern nucleus 6% of the 2.3 µm continuum maximum). For every panel, the flux is normalized in the zoomed FoV.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Line ratios diagram used to compute the extinction. Blue points are the observed line ratios of Paβ/Pαα and [Fe II]λ1.257 µm/[Fe II]λ1.644 µm. Dashed lines fit extinction values computed assuming theoretical line ratios based on the Cardelli et al. (1989) extinction law and different electronic temperatures and densities (as indicated). The symbols along the model lines mark optical extinction values in steps of $A_V = 1$, increasing from right to left.</figcaption></figure>

sensitive to line-emitting gas density and temperature, this is not true for the recombination lines (Ramos Almeida et al. 2009). The theoretical values of the two line ratios were taken from Hummer & Storey (1987) for Paschen lines and Bautista & Pradhan (1998) for [Fe II] lines. As shown in the figure, the model with $T_e = 3 \times 10^4$ K and $n_e = 10^5$ cm$^{-3}$ (i.e., red curve) reproduces best the data. In Fig. 8 we plot the extinction maps obtained from the Paschen decrement and the [Fe II] lines. As can be noticed, there are discrepancies between the reddening estimations derived with the two different line

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Extinction maps obtained from the iron line ratio (left panel) and the Paschen decrement (right panel) with the extinction law of Cardelli et al. (1989). Crosses indicate the positions of the nuclei.</figcaption></figure>

ratios, as discussed in Riffel et al. (2006): given that the intrinsic ratio of the recombination lines is likely influenced by high density and radiative transfer effects, it is difficult to determine which species more accurately reflects the actual extinction in the galaxy. Therefore, we used both line ratios to estimate the extinction. After observing that the dispersion from the red curve is due to the noise of the spaxels with the lowest S/N, we projected all the observed points perpendicularly to the model curve and calculated the extinction for every spaxel. We show the extinction map used in the following in Fig. 9.

We observed the maximum extinction of $A_V \approx 10$ near the southern nucleus. This indicates the highest concentration of dust, which extends from there to north and east directions. Regarding the northern nucleus, we found $A_V \approx 8$. We observed patchy dust located between the nuclei, which leads to a mean value of $A_V \approx 7$ in that region. Instead, in the more external regions, we find an $A_V$ range of 2–5. The fact that the higher values of extinction are measured at the locations of the AGN confirms that the variation of continuum peaks with wavelength is (at least in part) due to obscuration effects observed in Fig. 6. We point out the connection between the extinction gradient
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Extinction map used for dust correction of data. It was obtained as a combination of the extinction maps of the iron line ratio and the Paschen decrement of Fig. 8 (see text in Section 3.5 for details). Crosses indicate the positions of the nuclei. Black contours are the same as those in Fig. 6, and indicate the regions where most of the continuum emission originates.</figcaption></figure>

inside the continuum contours and the variation of continuum peaks with wavelength (observable in Fig. 6). Indeed, focusing on these regions, the less extincted parts correspond to the location of the peak in the blue wavelengths, while the more reddened zones coincide with the AGN.

Consistent extinction values have been obtained from the stellar continuum from pXPF fitting. We used spectra from 1 to 3 µm, masking out all the gas emission lines and fitting only the continuum without any additive/multiplicative polynomial, finding $A_V \approx 8$ on the nuclei position and $A_V \approx 6$ in the region between them. These extinction values are consistent with those provided by the gas.

Lutz et al. (2003) estimated the dust absorption toward NGC 6240 in the MIR-FIR spectral range from ISO SWS data with a FoV of 20'' × 20'', and found values of $A_V \approx 15-20$. The discrepancy with our results is probably due to the different spectral range used, rather than the different spatial regions. Our measurements only probe to the depth at which the nuclei become optically thick in the K band. If we are not probing all the way into the nuclei, measurements at longer wavelengths will result in higher inferred values of $A_V$. Using B-K colors from HST, Gerssen et al. (2004) estimated $A_V \ge 3-5$ at the position of the northern peak and the two southern peaks of the continuum at 1.2 µm in Fig. 6 (which they call N1, N2, and N3, respectively; see their Fig. 3). Another estimation is that by Pollack et al. (2007), who, with Keck II observation in K and H bands, measured extinction from stellar clusters scattered across the galaxy. They found $A_V = 8$ near the southern nucleus, and $A_V = 5$ at the location of the northern peak of the bluest continuum map in Fig. 6. Kollatschny et al. (2020) estimated extinction with the Balmer decrement and obtained values of $A_V \approx 10$ in the nuclear regions. All these estimates are consistent with ours, despite the different estimation methods used.

After the estimation of the extinction for every spaxel, we can now correct all observed fluxes for reddening. Going forward, we will use extinction-corrected flux values.

# 4. Results

In this section we present the results obtained from these NIR-Spec data, divided into moment maps, ISM properties and kinematics of gas. Through these sections, we identify some structures with distinct kinematical features. In particular, the present NIRSpec data show the presence of outflows associated with the nuclei, two rotating disks, and a complex molecular gas distribution. For sake of clarity in Fig. 10 (left panel) we sketch these structures that will be discussed in more detail in the following.

## 4.1. Moment maps

We show the moment maps of [Fe II]λ1.644 µm, Paα, and H₂ 1–0 S(1), obtained with our fit model, in Fig. 11. The zero moment maps – the left column of panels – show clear differences between the ionized and molecular gas components. This supports the notion that the molecular gas follows different kinematics than the ionized gas, as already hypothesized in our fitting procedure. Regarding the ionized lines, they have a similar distribution, even if [Fe II]λ1.644 µm seems to be less concentrated around the nuclei. The brightest line is Paα and its emission is mostly located at the nuclei, with the other lines being fainter by about one order of magnitude. However, for each species, the extended distribution of the gas is observed in the region between the nuclei, to the NW of the northern nucleus, and the SE and SW of the southern nucleus (corresponding to the bridge, bubble, Ext 1, and Ext 2 regions in Hermosa Muñoz et al. 2025, respectively). The velocity maps have values ranging between –300 km s⁻¹ and +300 km s⁻¹. We point out the presence of a region in all the maps with high redshifted velocities in the order of 200–300 km s⁻¹, in the center of the FoV (see Section 4.3.1 for a detailed discussion). The differences between ionized and molecular species are more evident looking at the most blueshifted velocities in the molecular gas map, that is dominant in the bottom half of the image. Regarding the velocity dispersion maps, the molecular hydrogen features the lowest σ, with a median value of 250 km s⁻¹ in the entire FoV, against the more turbulent other species with a median σ of 400 km s⁻¹. We also recognized part of the high-σ “V”-like shape structure observed by Hermosa Muñoz et al. (2025), which is starting from the SE of the southern nucleus and it extends toward the NE direction the NW one, passing through the southern nucleus (see their Fig. 6).

In Fig. 12 we plot the moment maps for the H₂ 1–0 S(1) emission line, but now separated into its narrow (top row) and broad (bottom row) components. In Appendix B, we plot the same moment maps for [Fe II]λ1.644 µm and Paα in Figs. B.1, B.2, respectively. The flux maps of the broad component of [Fe II]λ1.644 µm and Paα show notable features near both nuclei. While the broad emission component of both lines is brightest at the southern nucleus, both maps show an extended emission toward the S direction. Additionally, near the northern nucleus, there is an extended structure in the E direction. The respective velocity dispersion maps also show elevated values at these locations, especially near the southern nucleus where the dispersion reaches a maximum of ≈600 km s⁻¹, and forms a structure that extends along the N–S direction, coinciding with an ionization cone defined in Section 4.2.1 and depicted in Fig. 13. These observed properties, characterized by extended structures, high relative velocities, and significant velocity dispersions, can be naturally explained by the presence of outflows. As we will discuss in Section 4.2.1, the coincidence of the high
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Main kinematical features and gas distribution observed in the NIRSpec FoV (left panel), over-imposed on the flux map of Hα from NFM MUSE data (right panel) with schematic identification of arms and bubbles already known from the literature. The two ellipses indicate the rotating disk of stars, where the approaching region is in blue and the receding one is in red (Section 3.1). The molecular gas distribution is shown in light blue (Sections 4.1 and 4.3.1). The ionization cones of the two AGN (Sections 4.2.1 and 4.2.2) are indicated in orange. The light blue cone represents the molecular outflow near the southern nucleus (Section 4.3.1). Black crosses and the black ellipse indicate respectively the nuclei and the high-velocity molecular hydrogen blob positions.</figcaption></figure>

$\sigma$ structure with an ionization cone in the southern nucleus supports this interpretation, as observed in other local Seyfert galaxies (e.g., Mingozi et al. 2019).

The molecular hydrogen maps show a more complex kinematic behavior in the broad component. We point out two high velocity features: from the southern nucleus an approaching (blueshifted) outflow with $v \approx -600 \text{ km s}^{-1}$ in the SE direction, and from the northern nucleus a vertical component at $\sim 700 \text{ km s}^{-1}$, extending toward the S direction.

## 4.2. Interstellar medium properties

In this section we explore the properties of the ISM in the nuclear region of NGC 6240. We exploit the most used NIR diagnostics in the literature.

### 4.2.1. Near-infrared excitation diagnostic diagram

To understand the nature of the nuclear activity responsible for the excitation of the observed emission lines, we used the line ratios H<sub>2</sub> 1–0 S(1)/Bry and [Fe II]λ1.257 µm/Paβ (respectively H<sub>2</sub>/Bry and [Fe II]/Paβ, from now on), following Larkin et al. (1998), who observed a strong linear correlation in the log-log plot of the two line ratios. To distinguish among the possible mechanisms that regulate the intensities of these lines, we use the threshold values derived by Riffel et al. (2013), estimated using a large sample of galaxies previously classified:

- Star-forming galaxies (SFGs): [Fe II]/Paβ < 0.6 and H<sub>2</sub>/Bry < 0.4;

- AGN: 0.6 < [Fe II]/Paβ < 2 and 0.4 < H<sub>2</sub>/Bry < 6;

- Low ionization nuclear emission-line regions (LINERs) or shock: [Fe II]/Paβ > 2 and H<sub>2</sub>/Bry > 6.

We plot in the top panels of Fig. 13 the maps of these line ratios. As we can see in the bottom left panel, spaxels are divided into shock-LINER and AGN excitation dominated regions, while SFG excitation is not detected in any spaxel. To understand bet-

ter the distribution and the morphology of these regions, we applied the color-coding of the excitation diagram to plot an excitation map, Fig. 13 (bottom right panel). Most of the FoV is dominated by shock-LINER excitation, as can be expected by a dusty merging system such as NGC 6240, except for the areas around the two nuclei, which are clearly dominated by AGN excitation. In particular, we distinguish the two elongated structures that, starting around each nucleus, are extended in the S direction for the southern nucleus and in the NW for the northern one. There is also some AGN ionization in the east direction of the northern nucleus. We identify those structures with the ionization cones of the two AGN, which are sketched in the panel (see also Fig. 10). This result is in agreement with the excitation map of Riffel et al. (2021), who have a slightly wider FoV so they better appreciated and more clearly identified the northern structure (see their Fig. A1). Our H<sub>2</sub>/Bry map is also consistent with the H<sub>2</sub>/Paβ map by Medling et al. (2021) (see their Fig. 5).

Müller-Sánchez et al. (2011) and Medling et al. (2021) described the bubble and filamentary structures of NGC 6240 observed at a larger scale, in a FoV of almost 25'' × 25''. The most important ones are indicated in Fig. 10 (right panel). Müller-Sánchez et al. (2018) concluded that the NW bubble, detected in Hα with no [O III] detection, was possibly a starburst-driven outflow. The SE bubble is observed in both Hα and [O III], and they hypothesized a contribution of both AGN- and SF-ionization in this structure. Medling et al. (2021) found from BPT diagrams that all the arms and bubbles are in the LINER or shock-like ionization mechanism regime, except the E arm, which is more similar to the AGN-photo-ionization mechanism. Despite the differences in resolution and spatial scale between our data and Müller-Sánchez et al. (2011) and Medling et al. (2021) ones, the ionization cone of the northern nucleus we detected with the NIR diagnostic diagram (bottom right panel of Fig. 13) spatially aligns with the direction of the NW bubble. We are thus probably observing the launching region of this outflow, which would be AGN-driven in this
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Moment maps of [Fe II]λ1.644 µm, Paα, and H₂ 1–0 S(1) total emission lines with a threshold of S/N > 15. From left to right, we plot the flux map, the velocity map, and the velocity dispersion map. Crosses and the ellipse indicate, respectively, the nuclei and the high-velocity molecular hydrogen blob positions.</figcaption></figure>

scenario. This conclusion is also supported by the detection of coronal lines (CLs) in the same spatial direction (details in Section 4.2.2; see also Hermosa Muñoz et al. 2025). Similarly, the ionization cone observed near the southern nucleus and the molecular outflow in the SE direction (see Section 4.3.1) could be kinematically connected with the SE bubble. The connection of the kinematical structures discussed in this work and the ones observed on larger scale will be discussed by Tozzi et al. (in prep.) using a combination of IFU data from MUSE and KMOS on larger spatial scale. Regarding the other large-scale structures (i.e., the E arm and the W bubble), we did not find any correlation during the data analysis.

### 4.2.2. Detection of coronal line emission

Coronal lines result from collisionally excited forbidden transitions within low lying states of highly ionized species with ionization potentials IP > 100 eV and can be formed either in gas photoionized by a hard UV continuum or in a very hot collisionally ionized plasma. Indeed, they can be considered as indicative of AGN activity (e.g., Oliva et al. 1994, Moorwood et al. 1997). Most of these highly ionized lines have been observed in NIR spectra of nearby Seyfert galaxies (e.g., Rodríguez-Ardila et al. 2011, Lamperti et al. 2017, den Brok et al. 2022). While no comprehensive quantitative studies have been conducted on the

prevalence of CLs in the AGN population, some evidence suggests they are more likely observed in unobscured type-1 sources than in obscured ones. This aligns with recent detections of CLs in the nearby Seyfert 1.5 nucleus of NGC7469 using NIR-Spec and MIRI (U et al. 2022; Armus et al. 2023; Bianchin et al. 2024), contrasted with their absence in the highly obscured active nuclei of VV114 and Arp 220 (Perna et al. 2024).

Armus et al. (2006) detected the [Ne V]λ14.3 µm emission with Spitzer/IRS integrated spectra of NGC 6240. The other CL detection in NGC 6240 reported in literature is that of [Ca VIII]λ2.32 µm by Ilha et al. (2016). This CL is found at both nuclei as well as in the region between them, with the strongest emission observed at the southern nucleus.

In our NIRSpec data, we detected resolved emission of the following CLs only near the northern nucleus: [Si VI]λ1.96 µm (IP ~ 167 eV), [Si VII]λ2.48 µm (IP ~ 205 eV) and [Mg VIII]λ3.03 µm (IP ~ 225 eV). We fit these CLs with the same method presented in Section 3. Being Si lines blended with some H₂ lines, one of the two kinematically independent groups of the fitting procedure is the molecular family of H₂ 1–0 S(3), H₂ 1–0 Q(6) and H₂ 1–0 S(1)<fnref n="4" /> and the other one

<footnote n="4">⁴ We decided to fit also the H₂ 1–0 S(1) because it could help the fit to constrain the molecular line profile, being a strong and isolated emission line.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Moment maps of H₂ 1–0 S(1) narrow and broad components. The panels in the upper (lower) row correspond to the moment maps for the flux, velocity, and velocity dispersion of the narrow (broad) component. Crosses and the ellipse indicate respectively the nuclei and the high-velocity molecular hydrogen blob positions. We masked spaxels with S/N < 5.</figcaption></figure>

is the CLs group. Due to the low S/N of the CLs, we only use up to 2 Gaussian components. In Fig. 14 we plot the moment maps of the CLs in the northern nuclear region. The flux maps reveal that these species are predominantly concentrated in the nuclear region, with an extension toward the NW direction. The flux peak correspond to the northern nucleus position, confirming our assumption that Paα peaks at the AGN positions (see Section 2). All three CLs show elongated emission extending over ~0.7" (~350 pc) in the NW direction, which aligns closely with the ionizing cone contour identified in Section 4.2.1. Notably, [Si VI]λ1.96 µm, which has the lowest IP among the observed species, displays higher fluxes at greater distances from the nucleus. In contrast, [Mg VIII]λ3.03 µm, characterized by a higher IP, shows comparatively lower brightness. Furthermore, silicon lines are detected also in the detached structure of the ionization cone contour at almost 1.2" (~900 pc), in the top right angle of the FoV. This is consistent with the scenario that most energetic photons do not reach large distances from the AGN and are absorbed earlier.

The CLs velocity dispersion maps are similar to each other, with values ranging from 100 to 300 km s⁻¹. Qualitatively, the velocity maps are similar: the species present a blueshifted structure at 0.3" with ~−100–200 km s⁻¹ and a redshifted one at 0.6" with ~+300–400 km s⁻¹.

To detect these CLs in the southern nucleus – previously undetectable in individual spaxels around the nucleus – and to identify other faint CLs around both nuclei, we fit the spectra integrated over a radius of 5 spaxels on the nuclei position. This approach was used to enhance the signal-to-noise ratio (S/N). In addition to the previous CLs, we detected [Ca VIII]λ2.32 µm (IP ~ 127 eV) and [Si IX]λ2.58 µm (IP ~ 303 eV) in both nuclei

regions. Only [Si VI] was not detected in southern nucleus. In this case, we kinematically separate the CLs group from the other lines<fnref n="5" />. The fit spectra are shown in Fig. 15, where we plot the integrated spectra of the northern and southern nuclei. The CLs' fit luminosities are listed in Table 2.

We do not succeed in detecting a [Ca VIII]λ2.32 µm spatially resolved emission as done by Ilha et al. (2016) with Gemini NIFS data. In cube1 we do not detect any CLs. This is probably due to the higher absorption of dust in that datacube and the lower S/N of the spectra. The analysis of cube3 is beyond the scope of this work, but from the reduced datacube we can state the presence of other CLs in the northern nucleus, including [Ca IV]λ3.21 µm, [Al VI]λ3.66 µm, [Si IV]λ3.93 µm, [Mg IV]λ4.49 µm and [Ar VI]λ4.53 µm (see Figs. 1 and 2).

Hermosa Muñoz et al. (2025) studied NGC 6240 with JWST/MIRI data. Thanks to its spectral range of 5–28 µm and larger FoV, they detected some high excitation lines (such as [Mg V], [Ne V], and [Fe VIII]) in the N nucleus and in a bubble-like shape toward the NW direction of almost 5.2" (~2.7 kpc), outside of the NIRSpec FoV (see Fig. 6 in Hermosa Muñoz et al. 2025). This extended emission results completely redshifted and with a relatively low σ (~200 km s⁻¹), in accordance with our moment maps (see Fig. 14). This means that with NIRSpec we are looking at the launch region of this highly ionized bubble located out of our FoV. Regarding the southern nucleus, they detected [Mg V] and [Ne VI] in the integrated spectrum, after subtracting the strong PAH features.

<footnote n="5">⁵ We also fit Pfe because it's blended with [Mg VIII] in the integrated spectra, differently from the spatially resolved spectra.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Line ratio maps and diagnostic diagram for NIR emission lines. Top panels: [Fe II]λ1.257 µm/Paβ and H₂ 1–0 S(1)/Bry maps. Black contours delimit AGN-excited spaxels. Bottom-left panel: NIR diagnostic diagram. The dotted lines are the limit values dividing SF-, AGN-, and shock/LINER-excitation (see text for details). No SF-excitation points are detected. Bottom-right panel: Excitation map color coded according to the diagnostic diagram. Shock/LINER-excited spaxels are in green, and AGN-excited spaxels in red. Crosses and the ellipse indicate respectively the nuclei and the high-velocity molecular hydrogen blob positions. White contours are the same as those in Fig. 6, and indicate the regions where most of the continuum emission originates. Black cones indicate the approximate position of the two ionization cones. We used a S/N threshold of three, due to faint Bry emission.</figcaption></figure>

### 4.2.3. Excitation mechanisms and temperatures of H₂ lines

Following Mouri (1994), we investigated the dominant excitation mechanism of H₂ lines. Using H₂ line ratios it is possible to spatially separate the physical process undergoing in the galaxy. Molecular hydrogen IR emission lines can be due to two mechanisms:

- thermal process: collisions, when molecular clouds are heated by processes such as shocks, X-rays, or UV photons;

- non-thermal process: radiative decay resulting from the absorption of a UV photon in the Lyman-Werner band (912–1108 Å) or the collision with a fast electron due to X-ray ionization.

In each case, we expect a different emission-line spectrum, and therefore, the relative intensities among the H₂ emission lines may be used to discriminate between the powering mechanisms, keeping in mind that multiple emitting regions and multiple mechanisms may simultaneously be present.

Mouri (1994) suggested the use of the ratios 2–1 S(1) 2.247 µm/1–0 S(1) 2.121 µm versus 1–0 S(2) 2.033 µm/1–0 S(0) 2.223 µm. Another possible plot includes the line ratios 2–1 S(1) 2.247 µm/1–0 S(1) 2.121 µm versus 1–0 S(3) 1.958 µm/1–0 S(1) 2.121 µm: we did not use this latter because in our spectra the

1–0 S(3) line is blended with the [Si VI]λ1.963 µm, so its intensity cannot be reliable. In Fig. 16, we plot previous line ratios using the color-coding obtained in the previous section, which means that we separate the shock/LINER-dominated regions from the AGN-dominated ones. As one can expect, if we simply plot every spaxel we notice a high dispersion, assignable to spaxels with low S/N or those in the nuclear regions previously subjected to wiggles correction. Therefore, we performed a Voronoi binning to increase the S/N, using the faintest line, that is H₂ 1–0 S(0), to set a threshold at S/N = 4. We fit these new spectra with the same procedure described in Section 3 and we computed the H₂ line ratios to plot them in Fig. 16 as red and orange crosses, respectively for the southern and northern nuclei. We computed the median and standard deviation of all the points and crosses (black square with error bars). In the same panel, we plot a black curve representing the predictions for thermal emission, where the gas is in local thermodynamic equilibrium (LTE) and so the rotational and vibrational temperature have the same values (see next paragraphs for details). Each point on the curve represents temperature from 1000 to 6000 K, in steps of 1000 K, from left to right.

Fig. 16 also shows values of model calculations found in the literature, to better constrain the excitation mechanism.
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Moment maps around the northern nucleus of [Si VI]λ1.96 µm, [Si VII]λ2.48 µm, and [Mg VIII]λ3.03 µm emission lines. From left to right, we plot the flux map, the velocity map, and the velocity dispersion map. The cross indicates the northern nucleus position. The black contour represents the ionization cone of the northern nucleus from Fig. 13. We masked spaxels with S/N < 5.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Continuum-subtracted integrated spectra of the northern and southern nucleus in the upper and bottom panel, respectively. The black curve represents data, the red dotted curve represents the total model, the orange curve is the model of CLs, and the blue curve is the model of other lines.</figcaption></figure>

nisms involved. Black & van Dishoeck (1987) calculated UV-excited low-density gas models that occupied the orange ellipse on the right, while UV-excited high-density models by Sternberg & Dalgarno (1989) are represented by the green filled box on the left. The purple stars are models of X-ray-irradiated gas from Lepp & McCray (1983) and Draine & Woods (1990). Dors et al. (2012) calculated these ratios with the photoioniza-

tion code CLOUDY (Ferland 1993) for AGN models, pink polygon in the plot. The filled blue square is a UV thermal model by Davies et al. (2003) and the cyan star is from the shock model of Kwan et al. (1977).

As we can see in the figure, points in the diagram are away from the region predicted by the non-thermal UV excitation models of Black & van Dishoeck (1987) (identified by the
<table><caption>**Table 2.** High ionization lines detected in NGC 6240.</caption><thead><tr><th>Line</th><th>λ<sub>vac</sub><br>(Å)</th><th>IP<br>(eV)</th><th>log L<sub>line</sub><br>(erg s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>[Si VI]</td><td>19625</td><td>167</td><td>39.07 ± 0.05 (N)<br>not detected (S)</td></tr><tr><td>[Ca VIII]</td><td>23210</td><td>127</td><td>38.66 ± 0.06 (N)<br>38.8 ± 0.1 (S)</td></tr><tr><td>[Si VII]</td><td>24807</td><td>205</td><td>39.21 ± 0.04 (N)<br>39.42 ± 0.10 (S)</td></tr><tr><td>[Si IX]</td><td>25839</td><td>303</td><td>38.5 ± 0.8 (N)<br>38 ± 1 (S)</td></tr><tr><td>[Mg VIII]</td><td>30276</td><td>225</td><td>39.01 ± 0.13 (N)<br>39.32 ± 0.12 (S)</td></tr></tbody></table>

**Notes.** From left to right: Name of the CLs, rest-frame wavelengths, ionization potentials, and log luminosities in the northern (N) and southern (S) nucleus. See Section 4.2.2 for details.

<figure><img src="image_17.png" /></figure>

**Fig. 16.** H₂ excitation diagram. Green and red points are respectively the shock/LINER-dominated and the AGN-dominated spaxels (from Section 4.2.1). Red and orange crosses represent the rebinning regions near the southern and northern nuclei. Black square with error bars are the median and the standard deviation of all these points. Black circle with error bars are derived from the estimation of T<sub>rot</sub> and T<sub>vib</sub> by Riffel et al. (2021) integrating NGC 6240 data in a slightly wider FoV. Black dashed curve represents the thermal emission and on the secondary axis, we have T<sub>rot</sub> and T<sub>vib</sub>. Each point on the curve represents temperatures from 1000 to 6000 K, in steps of 1000 K, from left to right. The orange ellipse represents the region occupied by non-thermal UV excitation models of Black & van Dishoeck (1987) and the pink polygon covers the region of the photoionization code CLOUDY for AGN models of Dors et al. (2012). The green filled rectangle covers the locus of the thermal UV excitation models of Sternberg & Dalgarno (1989). The filled blue squares are the UV thermal models from Davies et al. (2003). The purple stars are models of X-ray-irradiated gas from Lepp & McCray (1983) and Draine & Woods (1990) and the cyan star is from the shock model of Kwan et al. (1977).

orange ellipse), indicating that the H₂ emission originates mainly in thermal processes. Indeed, our data are concentrated around the thermal emission at 2000 K but with a mild dispersion that could be an effect of low S/N. The models that better reproduce them are the ones of X-ray-irradiated gas and the AGN photoionization CLOUDY models. In particular, these latter represent the AGN-dominated regions – crosses and red points in Fig. 16–, as expected. We also observe that the ionization cones where

we rebinned, indicated by crosses in the plot, exhibit the most extreme values. This suggests that these regions are slightly farther from the thermal emission and may contain some contamination from non-thermal emission.

According to Rodríguez-Ardila et al. (2005), an alternative way to determine the mechanism driving the thermal or non-thermal origin of the H₂ excitation processes is to calculate the rotational temperature T<sub>rot</sub> and the vibrational temperature T<sub>vib</sub> of the molecular gas. Indeed, in a gas dominated by thermal excitation the two temperatures are quite similar, as expected for a gas in LTE. Non-thermal excitation, on the other hand, is characterized by a high vibrational temperature and a low rotational temperature; non-local UV photons, not characteristic of the local kinetic temperature, overproduce the highest energy levels compared to that expected for a Maxwell-Boltzmann population.

To estimate these temperatures we used the H₂ 2–1 S(1) 2.247 µm/H₂ 1–0 S(1) 2.121 µm and H₂ 1–0 S(2) 2.033 µm/H₂ 1–0 S(0) 2.223 µm line ratios, the same present in Fig. 16. We used N<sub>i</sub>/g<sub>i</sub> = F<sub>i</sub>λ<sub>i</sub>/(A<sub>i</sub>g<sub>i</sub>) = exp[-E<sub>i</sub>/(k<sub>B</sub>T<sub>exc</sub>)], where N<sub>i</sub> are the column densities in the upper level, g<sub>i</sub> the statistical weights, F<sub>i</sub> and λ<sub>i</sub> are the line fluxes and wavelengths, A<sub>i</sub> are the transition probabilities, E<sub>i</sub> are the energies of the upper level, T<sub>exc</sub> is the excitation temperature and k<sub>B</sub> is the Boltzmann constant. Using the transition probabilities from Turner et al. (1977), the rotational temperature is given by

$$T_{\text{rot}} \approx \frac{1113}{\ln\left(0.323 \frac{F_{\text{H}_2 \ 1-0 \text{S(2)}}}{F_{\text{H}_2 \ 1-0 \text{S(0)}}}\right)}, \quad (1)$$

and the vibrational temperature is given by

$$T_{\text{vib}} \approx \frac{5594}{\ln\left(1.355 \frac{F_{\text{H}_2 \ 1-0 \text{S(1)}}}{F_{\text{H}_2 \ 2-1 \text{S(1)}}}\right)}. \quad (2)$$

In Fig. 17 we show the maps of T<sub>rot</sub> and T<sub>vib</sub>, including the Voronoi regions. As noted in the previous paragraphs, the lowest T<sub>rot</sub> values and highest T<sub>vib</sub> values in the nuclear regions indicate non-thermal emission contamination near the nuclei. Comparing Fig. 16 with Fig. 17, the dispersion of the green points in the diagram is likely caused by extreme T<sub>rot</sub> values in spaxels near the FoV edges, which are probably affected by noise. The median values and standard deviations of these maps are T<sub>rot</sub> = (1894 ± 397) K and T<sub>vib</sub> = (2050 ± 141) K. Recent studies by Riffel et al. (2021), using data obtained with the NearInfrared Integral Field Spectrograph (NIFS; on the Gemini north Telescope) in the K band (slightly wider FoV and higher spectral resolution), computed T<sub>rot</sub> = (1247 ± 74) K and T<sub>vib</sub> = (2339 ± 3) K. Apart from the slight discrepancy with these latter results – possibly due to our higher sensitivity – we found a rotational and vibrational temperatures that are consistent with each other, as they should be in a gas in LTE. We conclude again, as done by using the diagnostic diagram of Fig. 16, that what we are seeing originates from thermal emission.

### 4.2.4. Mass of molecular gas

Following Riffel et al. (2008), we used the fluxes of the H₂ 1–0 S(1) and Paα emission lines to compute the mass of hot molecular and ionized hydrogen. The term “hot” depends on the fact that, compared to cold molecular gas (T < 100 K) detected via sub-millimeter CO line emission, IR H₂ roto-vibrational transitions trace a hotter phase of molecular gas with typical temperatures of T ~ 100–1000 K (e.g., García-Burillo et al. 2021).
<figure><img src="image_18.png" /><figcaption>Fig. 17. Logarithm of the rotational (left panel) and vibrational (right panel) temperatures computed from H₂ line ratios. Crosses indicate the positions of the nuclei.</figcaption></figure>

The mass of hot H₂ can be estimated, under the assumptions of local thermodynamic equilibrium and a T_vib = 2000 K, as proposed in Scoville et al. (1982), as

$$ \left(\frac{M_{\text{H}_2}}{M_{\odot}}\right) = 5.0776 \times 10^{13} \left(\frac{F_{\text{H}_2\text{1-0S(1)}}}{\text{erg s}^{-1} \text{cm}^{-2}}\right) \left(\frac{D}{\text{Mpc}}\right)^2, \quad (3) $$

where F<sub>H₂ 1-0S(1)</sub> is the dereddered H₂ 1-0 S(1) emission-line flux and D = 106 Mpc the galaxy distance (see Section 2). Integrating the flux over all the FoV, we obtained a value of M<sub>H₂</sub> ≈ 1.3 × 10<sup>5</sup> M<sub>☉</sub>. This value is in agreement with Riffel et al. (2021), where M<sub>H₂</sub> ≈ 1.1 × 10<sup>5</sup> M<sub>☉</sub>.

We compared the hot-to-cold ratios of molecular gas with the previous estimation by Medling et al. (2021). They found a ratio of ~10<sup>-6</sup> in the nuclear region between the nuclei, and roughly 100 times higher in the external filaments of the galaxy (see their Fig. 11). We also considered the results of Treister et al. (2020), who computed the cold molecular gas mass from CO (2-1) ALMA observations in three distinct regions indicated in Fig. 18: the sphere of influence of the northern and southern SMBHs and 1" diameter aperture centered equidistantly along the line connecting the two nuclei. They found 7.4 × 10<sup>8</sup> M<sub>☉</sub> for the northern nucleus, 3.3 × 10<sup>9</sup> M<sub>☉</sub> for the southern nucleus and 8.6 × 10<sup>9</sup> M<sub>☉</sub> for the central region. For a right comparison, we also calculated the hot molecular gas in those regions: we obtained 1.0 × 10<sup>4</sup> M<sub>☉</sub> for the northern nucleus, 2.3 × 10<sup>4</sup> M<sub>☉</sub> for the southern nucleus and 2.5 × 10<sup>4</sup> M<sub>☉</sub> for the central region. The respectively hot-to-cold ratios of molecular gas are therefore 1.4 × 10<sup>-5</sup>, 7.1 × 10<sup>-6</sup> and 2.9 × 10<sup>-6</sup>. These integrated values follow the same trend of the spatially resolved ratio of Medling et al. (2021), with lower values in the region between nuclei and higher values on the nuclei. Our estimations are well within the 10<sup>-7</sup>-10<sup>-5</sup> range found for molecular gas in starburst galaxies and AGN by Dale et al. (2005). In particular, the highest ratio of the northern nucleus is slightly lower but comparable to the extreme ratio of 6-7 × 10<sup>-5</sup> observed in the outflows of local LIRGs, for example by Emonts et al. (2014) and Pereira-Santaella et al. (2016).

### 4.2.5. Electron density

In order to compute the mass of ionized gas, we needed to estimate the electron density (see Section 4.2.6 and Eq. (6)). From a spectrum integrated into the same FoV of NIRSpec data, we fit the [S II]λ6716 Å and [S II]λ6731 Å emission lines observed in the MUSE NFM datacube<fnref n="6" />, with the same multi-

<figure><img src="image_19.png" /><figcaption>Fig. 18. Flux of the Paα emission line in logarithmic scale. Crosses indicate the positions of the nuclei. Circles indicate the spheres of influence of each SMBH and a 1" diameter aperture centered equidistantly along the line connecting the two nuclei used by Treister et al. (2020).</figcaption></figure>

Gaussian fitting procedure described in Section 3. We estimated a value of n<sub>e</sub> = (359 ± 52) cm<sup>-3</sup> using the module PyNeb (with T<sub>e</sub> = 10<sup>4</sup> K).

Another estimation of density can be achieved using some NIR [Fe II] lines, as described in Koo et al. (2016), such as [Fe II]λ1.534 µm, [Fe II]λ1.600 µm, [Fe II]λ1.644 µm, [Fe II]λ1.664 µm, and [Fe II]λ1.677 µm. Indeed, these are multiple lines of comparable excitation energies, with ratios that are mainly a function of electron density n<sub>e</sub> and depend weakly on temperature T<sub>e</sub>. We needed to significantly increase the S/N of the [Fe II] emission lines. In fact, except for the strong [Fe II]λ1.644 µm, other lines are almost undetected at spaxel level. We decide to extract an integrated spectrum from an aperture located in the FoV center with a radius of 30 spaxels – or 1.5"–, deriving [Fe II]λ1.677 µm and [Fe II]λ1.664 µm emission lines fluxes. After applying the usual fitting procedure, we computed the electron density with the equations

$$ \log n_e = 3.14 + 13.7 \frac{F_{1.664}}{F_{1.644}}, \quad \text{with} \quad 0.03 \le \frac{F_{1.664}}{F_{1.644}} \le 0.13; \quad (4) $$

$$ \log n_e = 2.93 + 8.14 \frac{F_{1.677}}{F_{1.644}}, \quad \text{with} \quad 0.07 \le \frac{F_{1.677}}{F_{1.644}} \le 0.23. \quad (5) $$

The two [Fe II] diagnostics gave comparable values, with a mean of log(n<sub>e</sub>/cm<sup>-3</sup>) = 4.12. This density exceeds that derived from optical [S II] lines, around 359 cm<sup>-3</sup>. A plausible interpretation of this is that the two methods trace different gas regions. This would mean that [Fe II] lines, with critical densities higher than [S II] (n<sub>c</sub> ~ 1600 cm<sup>-3</sup> for [S II]λ6716 Å; n<sub>c</sub> ~ 10<sup>4</sup> cm<sup>-3</sup> for iron lines), delineate post-shock regions experiencing greater compression. At the same time, due to different ionization energy (E<sub>ion</sub> = 7.9 eV for iron and E<sub>ion</sub> = 10.4 eV for sulfur), we expect that [S II] identifies regions with a higher degree of ionization. For comparison, we note that Hermosa Muñoz et al. (2025) found a density of n<sub>e</sub> ~ 1800 cm<sup>-3</sup>, using the ratio of [Ne V]λ14 µm and [Ne V]λ24 µm obtained from an integrated MIRI spectrum on the northern nucleus.

### 4.2.6. Mass of ionized gas

Similar to Riffel et al. (2021), we calculated the mass of ionized gas M<sub>HII</sub>, assuming an electron temperature T<sub>e</sub> = 10<sup>4</sup> K, with the

<footnote n="6">⁶ To reduce the degeneracy in the doublet, we fit this latter spectrum with Hα, Hβ, [S III]λ9069 Å lines and doublets of [O III]λλ5007,4959 Å, [N II]λλ6548,83 Å, and [O I]λλ6300,64 Å.</footnote>
<figure><img src="image_20.png" /><figcaption>Fig. 19. Logarithm ratio between the ionized gas mass and the hot molecular gas mass. The color bars depict the two values obtained with different proxy for the density estimation, namely [Fe II] and [S II] line ratios (see text for details). Crosses and the ellipse indicate respectively the nuclei and the high-velocity molecular hydrogen blob positions.</figcaption></figure>

formula

$$ \left(\frac{M_{HII}}{M_{\odot}}\right) = 2.30 \times 10^{18} \left(\frac{F_{Pao}}{\text{erg s}^{-1} \text{cm}^{-2}}\right) \left(\frac{D}{\text{Mpc}}\right)^2 \left(\frac{N_e}{\text{cm}^{-3}}\right)^{-1}. \quad (6) $$

where $F_{Pao}$ is the Pao integrated flux and $n_e$ is the electron density.

Using the value of $n_e \sim 359 \text{ cm}^{-3}$, obtained with the [S II] line ratio from MUSE NFM data (see Section 4.2.5), we obtained $M_{HII} \approx 2.8 \times 10^7 \text{ M}_{\odot}$, similar to the one obtained by Riffel et al. (2021) of $M_{HII} \approx 1.0 \times 10^7 \text{ M}_{\odot}$ (they assumed $n_e = 500 \text{ cm}^{-3}$). Using a higher density in Equation (6), we expect lower mass values. Indeed, using density from iron lines (see Section 4.2.5) we computed $M_{HII} \approx 7.7 \times 10^5 \text{ M}_{\odot}$. Given that we do not know which tracer is better suited to measure the density of the bulk of the gas, we will consider both extremes in the following.

We show in Fig. 19 the logarithm ratio between the ionized gas mass and the hot molecular gas mass. For the ionized mass, we used both the density estimations obtained above. The figure displays the ratio calculated with both iron lines and sulfur lines, as it is indicated by the color bars. We computed it considering all the components (narrow + broad) of the emission lines. Here it is clear the enhancement of ionized gas with respect to the molecular one in the two ionization cones, defined with the NIR diagnostic diagram in Section 4.2.1. Furthermore, the highest concentration of molecular gas with respect to the ionized one is located in the position of the highly redshifted blob of molecular hydrogen (see Section 4.3.1) indicated with a white ellipse. In the hypothesis that all the outflow is in the broad component, we calculated the same ratio $M_{H_2}/M_{HII}$ in the outflow present in our FoV. We found a logarithm value of 1.0 with iron lines density and 2.6 with sulfur lines density.

## 4.3. Gas kinematics

### 4.3.1. Molecular gas connecting the two nuclei

In Fig. 20 we plot channel maps of the emission line $H_2 1-0 \text{ S(1)}$, starting from about $-800 \text{ km s}^{-1}$ until $800 \text{ km s}^{-1}$. This figure contains a lot of interesting features: in the first row, in the bluer maps, an outflow can be recognized – highlighted with a white

contour in the figure – starting from the south nucleus, pointing in the SE direction. In the central row, a filamentary structure with three lobes can be observed, near the systemic velocity. From the right panel in the central row, the most interesting structure of the gas molecular kinematics appears increasingly in the channel maps [101 km s$^{-1}$, 265 km s$^{-1}$], where it seems there are two arms connecting the north nucleus and a region between the two nuclei. This region then continues to be connected to the north nucleus through an NS-directed arm and stays visible until about 800 km s$^{-1}$.

A significant reservoir of high-velocity molecular gas between the nuclei has been already seen in the cold molecular component (see Section 1): it was studied and discussed by Cicone et al. (2018) with [C I] (1–0) line and by Treister et al. (2020) with CO (2–1) transition, both with ALMA data. Also in Fyhrie et al. (2021) they observed this redshifted component, thanks to the channel maps of CO (3–2) and CO (6–5) obtained with ALMA (see their Figs. 4–5). To date, given its high velocity, it has been speculated by Treister et al. (2020) that this clumpy concentration of material is an outflow of molecular gas expelled from the nuclear region of the merger system. From the channel map of Fig. 20, this reservoir of molecular gas seems to be connected to both nuclei with some gas bridges. To better understand the nature of this clump and the kinematic structures of the molecular stream connecting the nuclei, we have used a position-velocity (p–v) diagram on the $H_2 1-0 \text{ S(1)}$ emission. As shown in the left panel of Fig. 21, the path of the diagram follows the direction of the previously mentioned southern outflow, and then in order it passes through the south nucleus, the reservoir of $H_2$ and the north nucleus. In the right panel of the Fig. 21, we display the p–v diagram, with vertical lines indicating the position of the landmarks along the path – namely, in increasing order of the offset, the southern nucleus, the high-velocity $H_2$ reservoir, and the northern nucleus. The molecular hydrogen is mainly concentrated on the southern nucleus, as already stated in other works (e.g., Ilha et al. 2016). Near the southern nucleus, we observe the molecular outflow component at a maximum velocity of about $-800 \text{ km s}^{-1}$. This feature corresponds to the structure with high blueshifted velocity and high velocity dispersion observed in the moment maps of the molecular broad component (see Fig. 12). In Fig. 10, we compare it with other observed kinematical structures, showing that the molecular outflow is not aligned with the ionization cone of the southern nucleus. This misalignment could have the same physical origin as the observed spatial anticorrelation between $H\alpha$ and CO(2–1) observed by Medling et al. (2021). They suggest that the hot ionized gas is flowing outwards along the path of least resistance (i.e., where denser molecular gas is not present). It is also interesting that our hot molecular outflow has the same direction of the CO flux map extension in the SE direction (see their Fig. 10). Such behavior is anyhow not surprising in a ULIRG, as it can also be observed, for example, in the merger galaxy Mrk 231, where the molecular and the neutral outflows have different axes (Rupke & Veilleux 2011; Feruglio et al. 2015). These misalignment between different outflow gas phases point toward a complex interplay of the two phenomena during the merger, or the presence of different outflow episodes, or a combination of the two. Another possibility is that the molecular outflow and the ionization cone are part of the same structure, with the molecular outflow forming a shell surrounding the ionized gas dominated cone but with only one side clearly visible, possibly due to a projection effect. Since we have no evidence of the other side of the molecular outflow, we consider only the detected outflow in the following discussion.

<footnote n="7">⁷ A movie showing the flux for every velocity channel is available online</footnote>
<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Channel maps of H₂ 1–0 S(1) emission lines. Velocity bins are indicated in the bottom left of every panel (in km s⁻¹). Crosses and the ellipse indicate respectively the nuclei and the high-velocity molecular hydrogen blob positions. White contour is the outflow launched from the southern nucleus.</figcaption></figure>

Another interesting feature is a gradient in velocity, which links the H₂ reservoir (which has a bulk velocity of ~500 km s⁻¹) and the two nuclei, with well-defined, steep behavior. Furthermore, between the H₂ reservoir and the northern nucleus, there is a fainter component, representing the arm visible in the central bottom panel of Fig. 20 linking the two structures. Although the highly redshifted molecular reservoir has been observed since Tecza et al. (2000), to date, nobody has ever clarified its nature. This unprecedented view with higher sensitivity and resolution using NIRSpec still does not allow us to reach a definite conclusion. However, we suggest two possible explanations.

First, it could be an outflow launched by the northern nucleus and oriented in the S direction. Indeed, the filaments in the [-117, 101] km s⁻¹ channel and the arms in the redshifted channels, point to what is expected in an expanding bubble. A similar structure was observed by Cresci et al. (2023) in the z = 1.59 prototypical obscured AGN XID2028, where, using the modeling outflows kinematics in AGN 3D (MOKA3D, see Marconcini et al. 2023) on the [O III] emission line, they assessed that the observed two emission-line arms and the cavity between them are attributable to an expanding bubble and a collimated outflow. This scenario has a critical point: from the

p–v diagram in Fig. 21 it is clear that this cloud is also connected to the southern nucleus, so we should conclude that the two outflows interact in the blob position. But in this case, we would expect an enhanced excitation from shock in that region traceable for example by [Fe II] emission, which is not observed (see Fig. 13). Ohyama et al. (2003) proposed a scenario where the kinematics of H₂ and the CO and H₂ intensity peaks might be explained by a collision between the off-nuclear molecular gas concentration and an external superwind outflow from the southern nucleus (see their Fig. 6). This is not consistent with our results, indeed they asserted a link between the off-nuclear molecular gas and the southern nucleus only, excluding any role of the northern nucleus.

Another plausible phenomenon is an inflow. Indeed, molecular gas inflows toward the central region in interacting and/or merging galaxies are expected, as observed in some ULIRGs by Herrera-Camus et al. (2020), to periodically replenish the molecular gas that is consumed by star formation or the SMBH, and/or ejected by outflows. With this idea, the high-redshift blob can be thought of as the result of the interaction between the two progenitor galaxies: during the merging, the molecular gas has been concentrated in the gravitational center of the system and
<figure><img src="image_22.png" /><figcaption>Fig. 21. Position-velocity diagram on the H₂ 1–0 S(1) emission. Left panel: Velocity channel map of H₂ 1–0 S(1) of ~400 km s⁻¹. In red we show the path used in the p–v diagram in the right panel. Crosses and the ellipse indicate respectively the nuclei and the high-velocity molecular hydrogen blob positions. Right panel: Position-velocity diagram of the H₂ 1–0 S(1) emission line. Vertical dashed lines indicate the position of the landmarks along the path (S is the southern nucleus, H₂ is the molecular blob, and N is the northern nucleus). The flux is in log scale.</figcaption></figure>

now is directly feeding the two AGN. The main problem of this scenario is understanding how this blob can have such a high velocity relative to the nuclei.

Without other information, the outflow and inflow are in principle indistinguishable. Thus, as already stated, we cannot draw any definite conclusions.

### 4.3.2. Outflow energetics of molecular and ionized gas

Mass-outflow rates are often calculated based on physical parameters such as outflow mass ($M_{out}$), radius ($R_{out}$), and velocity ($v_{out}$). This calculation typically assumes a spherical or (multi-)cone geometry with a constant velocity throughout the volume, under the simplifying assumption that there is no temporal or spatial acceleration or deceleration of the flow. While this assumption may oversimplify the dynamics, neglecting potential effects such as ballistic deceleration, radiation pressure acceleration, or interaction with other winds, it represents a simplified but useful framework for interpreting current observations. Therefore, even though gas kinematics is clearly much more complex in the nuclear region of NGC 6240, we made these simplifying assumptions using a radius of almost half of the FoV of our observations. The outflow rate is then derived via

$$ \dot{M}_{\text{out}} = C \frac{M_{\text{out}} v_{\text{out}}}{R_{\text{out}}}. \qquad (7) $$

The value of C depends on the adopted outflow history (Lutz et al. 2020): assuming that the outflow has started at a point in the past at $-t_{\text{flow}} = -R_{\text{out}}/v_{\text{out}}$ and has continued with a constant mass outflow rate until now directly leads to $C = 1$. We computed the average outflow mass rate, the outflow energy, and the outflow power as

$$ \begin{aligned} \dot{M}_{\text{out}} &= \frac{M_{\text{out}} v_{\text{out}}}{R_{\text{out}}}, \\ E_{\text{out}} &= \frac{1}{2} M_{\text{out}} v_{\text{out}}^2, \\ \dot{E}_{\text{out}} &= \frac{1}{2} \dot{M}_{\text{out}} v_{\text{out}}^2. \end{aligned} \qquad (8) $$

We applied Equations (8) in a circular region with radius $R_{out}$ of 30 spaxels, that is 1.50", which represents almost all NIRSpect

FoV – excluding bad spaxels at the edge of the image. We calculated the $M_{out}$ for the broad component of Pa$\alpha$ (with both sulfur lines density and iron lines density) and H₂ 1–0 S(1) emission line in the following way: we used the broad component flux of the two species, described in Section 4.1, and we computed the masses with Equations (6) and (3) (see Sections 4.2.6 and 4.2.4). In the literature, there is not a single definition of the bulk velocity of the outflow (see the review by Harrison et al. 2018). We computed the outflow velocity as

$$ v_{\text{out}} = \max (|v_{10}^{\text{max}} - v_{\text{sys}}|, |v_{90}^{\text{max}} - v_{\text{sys}}|), \qquad (9) $$

where $v_{sys}$ is the systemic velocity of the galaxy (estimated from the redshift in Section 2 and set to 0 km s⁻¹) and $v_{10}^{max}$ and $v_{90}^{max}$ are the maximum value, respectively, of 10th and 90th velocity percentiles of the overall emission line profile (i.e., narrow + broad components). The necessity of this definition, widely adopted in the literature (e.g., Cano-Díaz et al. 2012; Cresci et al. 2015; Carniani et al. 2015; Marasco et al. 2020; Tozzi et al. 2021, 2024), arises from the unknown geometry and orientation of the outflow relative to the line of sight. Due to our lack of knowledge regarding the true angle of the outflow with respect to the observer’s line of sight, and considering that the majority of the outflow is unlikely to be directed toward the observer, we assumed that the most accurate representation of the outflow speed is provided by the velocity “tail” of the line profile, that is $v_{10}$ and $v_{90}$ as described in Eq. (9). As discussed in other studies (e.g., Cresci et al. 2015, Tozzi et al. 2024), these values are considered more appropriate for representing the outflow velocity than the mean or median velocity of the line, which can be significantly influenced by projection effects and dust absorption, this latter being very strong in NGC 6240. Thus we used $v_{out} \approx 915$ km s⁻¹ for Pa$\alpha$ emission and $v_{out} \approx 800–900$ km s⁻¹ for H₂ 1–0 S(1) emission. The results are presented in Table 3. For the ionized phase, we found a mass outflow rate in the range 0.2–9 M$_{\odot}$ yr⁻¹, depending on the assumption on the electron density. For the hot molecular gas, we found outflow rate of ~0.02 M$_{\odot}$ yr⁻¹.

These values are mostly indicative because we did not recognize a classical outflow conical shape across all FoV, but we calculated the properties of all the broad components in our
<table><caption>**Table 3.** Properties of the ionized and molecular outflow in NGC 6240.</caption><thead><tr><th>Species and spatial region</th><th>v<sub>out</sub> (km s<sup>-1</sup>)</th><th>M&#x0307;<sub>out</sub> (M<sub>&#x210D;</sub> yr<sup>-1</sup>)</th><th>E<sub>out</sub> (erg)</th><th>E&#x0307;<sub>out</sub> (erg s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Paα (with n<sub>e</sub><sup>[S II]</sup>), circular region</td><td>914.5</td><td>9.4</td><td>6.3 &times; 10<sup>55</sup></td><td>2.5 &times; 10<sup>42</sup></td></tr><tr><td>Paα (with n<sub>e</sub><sup>[Fe II]</sup>), circular region</td><td>914.5</td><td>0.26</td><td>1.7 &times; 10<sup>54</sup></td><td>6.8 &times; 10<sup>40</sup></td></tr><tr><td>H<sub>2</sub> 1–0 S(1), circular region</td><td>796.5</td><td>0.020</td><td>1.1 &times; 10<sup>53</sup></td><td>3.9 &times; 10<sup>39</sup></td></tr><tr><td>H<sub>2</sub> 1–0 S(1), SE outflow</td><td>787.3</td><td>0.003</td><td>1.4 &times; 10<sup>52</sup></td><td>6.4 &times; 10<sup>38</sup></td></tr></tbody></table>

**Notes.** Velocity, mass rate, kinetic energy, and power of the outflow for ionized and molecular gas. These values were calculated for a circular region centered in our FoV. In addition, we calculated the outflow energetics for the outflow near the southern nucleus in the SE direction, discussed in Section 4.3.1.

FoV, associating it with a pure outflow. As stated in the previous section in Fig. 20, in the blue channel maps of H₂ 1–0 S(1) emission line we could outline a collimated outflow in the SE direction starting from the southern nucleus. So we calculated the outflow energetics also for this structure, using $R_{\text{out}} = 550$ pc and $v_{\text{out}} \approx 790$ km s⁻¹, with the resulting values in Table 3. Cicone et al. (2018) estimated the outflow mass rate of the cold molecular component on a scale of 6 kpc, using the flux of CO (1–0) and obtaining a value of almost 2500 M☉ yr⁻¹. Despite the extension of the outflow, they stated that the largest contribution is given by the central region. Indeed, Treister et al. (2020) – with a FoV of about 6′′ × 6′′, more similar to ours – also found a compatible value with CO (2–1) ALMA observations. It is interesting to note that if we use our most conservative estimation of hot-to-cold ratio (assuming this ratio is similar for the broad component; see Section 4.2.4) for molecular gas of ~10⁻⁵, the estimation of the outflow mass rate for the hot molecular gas of 0.020 M☉ yr⁻¹ is in agreement with the cited literature's values for the cold component. Furthermore, using the smaller value of ~10⁻⁶, we found that the outflow mass rate for the cold molecular component should reach ~10⁴ M☉ yr⁻¹. Regarding the ionized component, its outflow rate is greater than the hot molecular one by 1–3 orders of magnitude, depending on the used value of electron density (see Section 4.2.5). However, we confirmed that the total molecular phase accounts for the majority of the outflowing gas once accounting for the cold phase, featuring masses approximately two orders of magnitude larger than those observed in ionized outflows, as observed by Fiore et al. (2017).

Fiore et al. (2017) identified a strong correlation between mass outflow rate, outflow power, and bolometric luminosity in AGN. As shown in Fig. 22, they reported $\dot{M} \propto L_{\text{bol}}^{0.76 \pm 0.06}$ and $\dot{E}_{\text{out}} \propto L_{\text{bol}}^{1.29 \pm 0.08}$ for molecular winds (blue dashed lines in the top and bottom panel, respectively), while for ionized winds the scaling is $\dot{M} \propto L_{\text{bol}}^{1.29 \pm 0.38}$ and $\dot{E}_{\text{out}} \propto L_{\text{bol}}^{1.48 \pm 0.37}$ (green dashed lines in top and bottom panel, respectively). Yamada et al. (2021) found an X-ray luminosity of $1.37 \times 10^{10}$ L☉ and $5.21 \times 10^{9}$ L☉ for the southern and northern nucleus, respectively. Applying the bolometric correction from Duras et al. (2020), we obtained $L_{\text{bol}} = 2.38 \times 10^{11}$ L☉ for the southern nucleus and $L_{\text{bol}} = 8.53 \times 10^{10}$ L☉ for the northern one. To calculate the total bolometric luminosity, we summed the values of the two nuclei. Using these values, along with those listed in Table 3⁸, we plot in Fig. 22 the outflow properties of the cold molecular gas within the circular region (black crosses) and the SE outflow (purple crosses), as well as the ionized gas within the circular region, employing both density estimations (orange crosses for the [Fe II] method

and red crosses for the [S II] method). For the calculation of the cold molecular gas mass, we assumed a conservative hot-to-cold ratio of 10⁻⁵. The ionized outflow values are consistent with both relations. Regarding the cold molecular ones, it seems that the values of the circular region are too high. This discrepancy is likely due to the inclusion of features in the broad component of H₂ that are associated with the merger rather than the outflow. This could suggest that the H₂ blob (discussed in Section 4.3.1) may not be an outflow but rather a tidal feature. Another possibility is that the hot-to-cold molecular gas ratio derived for the bulk of the gas is too low for the outflow phase (e.g., see Ulivi et al. 2025, who found in Arp 220 hot-to-cold molecular gas ratios of $2 \times 10^{-4}$ and $9 \times 10^{-5}$ in the outflows of the east and west nuclei, respectively), where a larger part of the cold gas could be destroyed or is undetectable in the NIR transitions probed by NIRSpec.

## 5. Conclusion

In this work, we have analyzed high spectral resolution data obtained from JWST/NIRSpec of the ULIRG NGC 6240, a local prototypical merger known to host at least two AGN. Despite the high sensitivity and spatial resolution of our NIRSpec observations, we did not find any evidence of the presence of the third nucleus proposed by Kollatschny et al. (2020). In Fig. 10, we show a graphical summary of the more interesting features observed with this data. In the following, we summarize our main findings.

We first subtracted the stellar continuum and then conducted an emission line multi-component fitting, dividing the emission species into two families, namely molecular and ionized lines, under the hypothesis that the two gas phases follow different kinematics. At the end of the fitting procedure, we also separated the broad and the narrow components in order to distinguish between different kinematic structures in the galaxy. Then we performed the dust extinction correction on the datacubes.

As for the physical properties of the ISM, our main results are the following:

- We studied the properties of the ISM using NIR emission-line diagnostics. From the H₂ 1–0 S(1)/Bry and [Fe II]λ1.257 µm/Paβ line ratios, we concluded that the circumnuclear regions of this galaxy are mostly excited from shock-LINER excitation, apart from the AGN excitation that delineates two ionization cones starting from the two nuclei (orange cones in Fig. 10).

- We detected spatially resolved emission of CLs around the northern nucleus. The emission peaks on the nucleus position and is extended in the NW direction, aligning closely with the ionizing cone found with NIR diagnostics. Using

<footnote n="8">⁸ Fiore et al. (2017) calculated outflow properties with $C = 3$ (see Eq. (7)), so we applied this factor to the values in the table.</footnote>
<figure><img src="image_23.png" /><figcaption>**Fig. 22.** Outflow energetics relations compared with the literature. Top panel: Mass outflow rate as a function of the AGN bolometric luminosity. Green points and green dashed line are the outflow values and correlation for the ionized gas from Fiore et al. (2017), while blue points and blue dashed line are those for the molecular gas. Black cross is the cold molecular outflow in the circular region and purple cross that in the SE outflow. Orange and red crosses are the ionized outflow in the circular region computed with the [Fe II] method and the [S II] method, respectively. Bottom panel: Outflow power as a function of the AGN bolometric luminosity. Solid, dashed and dotted line represent the correlations $\dot{E}_{\text{out}} = 1, 0.1, 0.01 L_{\text{bol}}$.</figcaption></figure>

integrated spectra, we detected CLs also in the southern nucleus.

- We investigated the excitation mechanisms of the molecular hydrogen lines using specific $H_2$ line ratios. We found that the nuclear region is mainly characterized by thermal emission at ~2000 K, with a non-thermal contribution near the two nuclei. We confirmed this result by computing the vibrational and rotational temperatures in our FoV, which are equal to $T_{\text{rot}} = (1894 \pm 397)$ K and $T_{\text{vib}} = (2050 \pm 141)$ K.

- We computed the mass of hot $H_2$ in the FoV, assuming local thermodynamic equilibrium and a $T_{\text{vib}} = 2000$ K and obtaining a mass of ~1.3 × $10^5 M_{\odot}$. Using previous literature results, we estimated the hot-to-cold ratio in some spatial regions of our FoV, finding values of $10^{-6}–10^{-5}$, which are within the range found in the literature. Regarding the ionized gas mass, the estimation depends on the electron density. For this, we adopted the value obtained with the [S II] line ratio and that with the [Fe II]. We obtained $M_{\text{HII}} \approx 2.8 \times 10^7 M_{\odot}$ and $M_{\text{HII}} \approx 7.7 \times 10^5 M_{\odot}$, respectively.

Regarding the kinematics of the stellar and ISM components our findings are as follows:

- We studied in detail the kinematics of the ISM and in particular that of the ionized and hot molecular gas components. We

presented the moment maps of the brightest emission lines for both the total component and the separated narrow and broad components. We highlighted the different kinematics of the molecular hydrogen with respect to the other lines, in particular its higher velocity values and different spatial distribution (light blue region in Fig. 10). We also found a match around the southern nucleus between the velocity map of the stars and that of the narrow component of $Pa\alpha$, meaning that the narrow component of this line follows the rotating disk motion in that region.

- From the channel maps, we discovered features never observed before with this resolution. A blue outflow at high velocity, ~800 km s⁻¹, is observable in the SE direction near the southern nucleus (light blue outflow in Fig. 10). Around the systemic velocity, there are some filaments that, with increasing velocity until ~400 km s⁻¹, assume a more complex structure, with some arms connecting the northern nucleus and a blob located between the nuclei. The most peculiar characteristic is that this blob remains observable up to ~800 km s⁻¹. Although its true nature remains ambiguous, we considered both the possibility of an outflow and an inflow.

- In the end, we estimated the energetics of the outflow in our FoV. We calculated the outflow mass rate, its kinetic energy, and its power for both the ionized gas and the molecular gas, finding that the molecular outflow largely dominates the outflow rate. In addition, we estimated these properties for the blue outflow observed near the southern nucleus.

In the future, we will analyze cube3, the less extinction-affected datacube, to better constrain the physical conditions in NGC 6240 by exploiting, for instance, PAH and CL emission. Observations in other bands with NIRSpec-like spatial and spectral resolution might allow us to fully understand the nature of the $H_2$ cloud with a parallel theoretical approach where simulations could explain its current kinematics. To conclude, these NIRSpec IFU observations have allowed us to characterize with unprecedented detail a complex system such as NGC 6240 and have demonstrated the fundamental role of JWST in the study of local nearby galaxies.

## Data availability

A movie is available at https://www.aanda.org

**Acknowledgements.** We thank the anonymous referees for their feedback and suggestions, which have helped improve this manuscript. GC, AM, MP and LU acknowledge the support of the INAF Large Grant 2022 “The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities”. GC and AM acknowledge support from PRIN-MUR project “PROMETEUS” (202223XPZM) and INAF Large Grant 2022 “Dual and binary SMBH in the multi-messenger era”. SA and MP acknowledge support from the research project PID2021-127718NB-I00 of the Spanish Ministry of Science and Innovation/State Agency of Research (MCIN/AEI/10.13039/501100011033). MP acknowledges support through the grant RYC2023-044853-I, funded by MICIU/AEI/10.13039/501100011033 and FSE+. KF acknowledges support through the ESA research fellowship programme. IL acknowledges support from PID2022-140483NB-C22 funded by AEI 10.13039/501100011033 and BDC 20221289 funded by MCIN by the Recovery, Transformation and Resilience Plan from the Spanish State, and by NextGenerationEU from the European Union through the Recovery and Resilience Facility, and from PRIN-MUR project “PROMETEUS” financed by the European Union – Next Generation EU, Mission 4 Component 1 CUP B53D23004750006. AJB acknowledges funding from the “FirstGalaxies” Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 789056). GT acknowledges financial support from the European Research Council (ERC) Advanced Grant under the European Union’s Horizon Europe research and innovation programme (grant agreement AdG GALPHYS, No. 101055023).
References

Alexander, D. M., & Hickox, R. C. 2012, New A Rev., 56, 93

Armus, L., Bernard-Salas, J., Spoon, H. W. W., et al. 2006, ApJ, 640, 204

Armus, L., Lai, T., U, V., et al. 2023, ApJ, 942, L37

Baan, W. A., Hagiwara, Y., & Hofner, P. 2007, ApJ, 661, 173

Bautista, M. A., & Pradhan, A. K. 1998, ApJ, 492, 650

Bennett, C. L., Larson, D., Weiland, J. L., & Hinshaw, G. 2014, ApJ, 794, 135

Bianchin, M., U, V., Song, Y., et al. 2024, ApJ, 965, 103

Black, J. H., & van Dishoeck, E. F. 1987, ApJ, 322, 412

Böker, T., Arribas, S., Lützgendorf, N., et al. 2022, A&A, 661, A82

Böker, T., Beck, T. L., Birkmann, S. M., et al. 2023, PASP, 135, 038001

Cano-Díaz, M., Maiolino, R., Marconi, A., et al. 2012, A&A, 537, L8

Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345

Cappellari, M., & Emsellem, E. 2004, PASP, 116, 138

Cappellari, M., Emsellem, E., Bacon, R., et al. 2007, MNRAS, 379, 418

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Carniani, S., Marconi, A., Maiolino, R., et al. 2015, A&A, 580, A102

Cicone, C., Severgnini, P., Papadopoulos, P. P., et al. 2018, ApJ, 863, 143

Cresci, G., Marconi, A., Zibetti, S., et al. 2015, A&A, 582, A63

Cresci, G., Tozzi, G., Perna, M., et al. 2023, A&A, 672, A128

Dale, D. A., Sheth, K., Helou, G., Regan, M. W., & Hüttemeister, S. 2005, AJ, 129, 2197

Davies, R. I., Sternberg, A., Lehnert, M., & Tacconi-Garman, L. E. 2003, ApJ, 597, 907

den Brok, J. S., Koss, M. J., Trakhtenbrot, B., et al. 2022, ApJS, 261, 7

D'Eugenio, F., Pérez-González, P. G., Maiolino, R., et al. 2024, Nat. Astron., 8, 1443

Dors, O. L., Jr., Riffel, R. A., Cardaci, M. V., et al. 2012, MNRAS, 422, 252

Draine, B. T., & Woods, D. T. 1990, ApJ, 363, 464

Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73

Emonts, B. H. C., Piqueras-López, J., Colina, L., et al. 2014, A&A, 572, A40

Engel, H., Davies, R. I., Genzel, R., et al. 2010, A&A, 524, A56

Fabbiano, G., Paggi, A., Karovska, M., et al. 2020, ApJ, 902, 49

Fahrion, K., Böker, T., Perna, M., et al. 2024, A&A, 687, A83

Ferland, G. J. 1993, Hazy, A Brief Introduction to Cloudy 84, 565

Feruglio, C., Fiore, F., Carniani, S., et al. 2015, A&A, 583, A99

Fiore, F., Feruglio, C., Shankar, F., et al. 2017, A&A, 601, A143

Fyhrie, A., Glenn, J., Rangwala, N., et al. 2021, ApJ, 922, 208

Gallimore, J. F., & Beswick, R. 2004, AJ, 127, 239

García-Burillo, S., Alonso-Herrero, A., Ramos Almeida, C., et al. 2021, A&A, 652, A98

Genzel, R., Lutz, D., Sturm, E., et al. 1998, ApJ, 498, 579

Gerssen, J., van der Marel, R. P., Axon, D., et al. 2004, AJ, 127, 75

Harrison, C. M., Costa, T., Tadhunter, C. N., et al. 2018, Nat. Astron., 2, 198

Heckman, T. M., Armus, L., & Miley, G. K. 1990, ApJS, 74, 833

Hermosa Muñoz, L., Alonso-Herrero, A., Labiano, A., et al. 2025, A&A, 693, A321

Herrera-Camus, R., Sturm, E., Graciá-Carpio, J., et al. 2020, A&A, 633, L4

Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1

Hummer, D. G., & Storey, P. J. 1987, MNRAS, 224, 801

Ilha, G. d. S., Bianchin, M., & Riffel, R. A. 2016, Ap&SS, 361, 178

Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80

Joseph, R. D., Wade, R., & Wright, G. S. 1984, Nature, 311, 132

Kollatschny, W., Weilbacher, P. M., Ochmann, M. W., et al. 2020, A&A, 633, A79

Komossa, S., Burwitz, V., Hasinger, G., et al. 2003, ApJ, 582, L15

Koo, B.-C., Raymond, J. C., & Kim, H.-J. 2016, J. Korean Astron. Soc., 49, 109

Krajnović, D., Emsellem, E., Cappellari, M., et al. 2011, MNRAS, 414, 2923

Kwan, J. H., Gatley, I., Merrill, K. M., Probst, R., & Weintraub, D. A. 1977, ApJ, 216, 713

Lamperti, I., Koss, M., Trakhtenbrot, B., et al. 2017, MNRAS, 467, 540

Larkin, J. E., Armus, L., Knop, R. A., Soifer, B. T., & Matthews, K. 1998, ApJS, 114, 59

Law, D. R., E. Morrison, J., Argyriou, I., et al. 2023, AJ, 166, 45

Lepp, S., & McCray, R. 1983, ApJ, 269, 560

Lípari, S., Terlevich, R., Díaz, R. J., et al. 2003, MNRAS, 340, 289

Lonsdale, C. J., Farrah, D., & Smith, H. E. 2006, in Astrophysics Update 2, ed. J. W. Mason, 285

Lutz, D., Sturm, E., Genzel, R., et al. 2003, A&A, 409, 867

Lutz, D., Sturm, E., Janssen, A., et al. 2020, A&A, 633, A134

Marasco, A., Cresci, G., Nardini, E., et al. 2020, A&A, 644, A15

Maraston, C., & Strömbäck, G. 2011, MNRAS, 418, 2785

Marconcini, C., Marconi, A., Cresci, G., et al. 2023, A&A, 677, A58

Max, C. E., Canalizo, G., Macintosh, B. A., et al. 2005, ApJ, 621, 738

Max, C. E., Canalizo, G., & de Vries, W. H. 2007, Science, 316, 1877

Medling, A. M., Kewley, L. J., Calzetti, D., et al. 2021, ApJ, 923, 160

Mingozzi, M., Cresci, G., Venturi, G., et al. 2019, A&A, 622, A146

Moorwood, A. F. M., Marconi, A., van der Werf, P. P., & Oliva, E. 1997, Ap&SS, 248, 113

Mouri, H. 1994, ApJ, 427, 777

Müller-Sánchez, F., Prieto, M. A., Hicks, E. K. S., et al. 2011, ApJ, 739, 69

Müller-Sánchez, F., Nevin, R., Comerford, J. M., et al. 2018, Nature, 556, 345

Nardini, E., Risaliti, G., Salvati, M., et al. 2009, MNRAS, 399, 1373

Nardini, E., Wang, J., Fabbiano, G., et al. 2013, ApJ, 765, 141

Ohyama, Y., Yoshida, M., & Takata, T. 2003, AJ, 126, 2291

Oliva, E., Salvati, M., Moorwood, A. F. M., & Marconi, A. 1994, A&A, 288, 457

Pereira-Santaella, M., Colina, L., García-Burillo, S., et al. 2016, A&A, 594, A81

Perna, M., Arribas, S., Marshall, M., et al. 2023, A&A, 679, A89

Perna, M., Arribas, S., Lamperti, I., et al. 2024, A&A, 690, A171

Pollack, L. K., Max, C. E., & Schneider, G. 2007, ApJ, 660, 288

Ramos Almeida, C., Pérez García, A. M., & Acosta-Pulido, J. A. 2009, ApJ, 694, 1379

Riffel, R., Rodríguez-Ardila, A., & Pastoriza, M. G. 2006, A&A, 457, 61

Riffel, R. A., Storchi-Bergmann, T., Winge, C., et al. 2008, MNRAS, 385, 1129

Riffel, R., Rodríguez-Ardila, A., Aleman, I., et al. 2013, MNRAS, 430, 2002

Riffel, R. A., Storchi-Bergmann, T., Riffel, R., et al. 2021, MNRAS, 504, 3265

Risaliti, G., Sani, E., Maiolino, R., et al. 2006, ApJ, 637, L17

Rodríguez-Ardila, A., Riffel, R., & Pastoriza, M. G. 2005, MNRAS, 364, 1041

Rodríguez-Ardila, A., Prieto, M. A., Portilla, J. G., & Tejeiro, J. M. 2011, ApJ, 743, 100

Rupke, D. S. N., & Veilleux, S. 2011, ApJ, 729, L27

Sanders, D. B., & Mirabel, I. F. 1996, ARA&A, 34, 749

Sanders, D. B., Soifer, B. T., Elias, J. H., Neugebauer, G., & Matthews, K. 1988, ApJ, 328, L35

Scoville, N. Z., Hall, D. N. B., Ridgway, S. T., & Kleinmann, S. G. 1982, ApJ, 253, 136

Scoville, N. Z., Evans, A. S., Dinshaw, N., et al. 1998, ApJ, 492, L107

Scoville, N., Sheth, K., Walter, F., et al. 2015, ApJ, 800, 70

Silk, J. 2013, ApJ, 772, 112

Souchay, J., Mathis, S., & Tokieda, T. 2013, Tides in Astronomy and Astrophysics (Berlin, Heidelberg: Springer-Verlag)

Sternberg, A., & Dalgarno, A. 1989, ApJ, 338, 197

Tacconi, L. J., Genzel, R., Tecza, M., et al. 1999, ApJ, 524, 732

Tecza, M., Genzel, R., Tacconi, L. J., et al. 2000, ApJ, 537, 178

Tozzi, G., Cresci, G., Marasco, A., et al. 2021, A&A, 648, A99

Tozzi, G., Cresci, G., Perna, M., et al. 2024, A&A, 690, A141

Treister, E., Messias, H., Privon, G. C., et al. 2020, ApJ, 890, 149

Turner, J., Kirby-Docken, K., & Dalgarno, A. 1977, ApJS, 35, 281

U, V., Sanders, D. B., Mazzarella, J. M., et al. 2012, ApJS, 203, 9

U, V., Lai, T., Bianchin, M., et al. 2022, ApJ, 940, L5

Ulivi, L., Perna, M., Lamperti, I., et al. 2025, A&A, 693, A36

van der Werf, P. P., Genzel, R., Krabbe, A., et al. 1993, ApJ, 405, 522

Veilleux, S., Shopbell, P. L., Rupke, D. S., Bland-Hawthorn, J., & Cecil, G. 2003, AJ, 126, 2185

Veilleux, S., Meléndez, M., Sturm, E., et al. 2013, ApJ, 776, 27

Wang, J., Nardini, E., Fabbiano, G., et al. 2014, ApJ, 781, 55

Yamada, S., Ueda, Y., Tanimoto, A., et al. 2021, ApJS, 257, 61
## Appendix A: Cube1 stellar subtraction

In this appendix, we present the analog examples of Fig. 3 for cube1. In Fig. A.1 we plot the fit cube1 spectra from the same spaxel of Fig. 3, one centered on the southern nucleus (upper panel) and the other located in the region between the two nuclei (lower panel). In both, we show data (in black), the total model of emission lines and stellar component (in red), the model of only stars (in green), and the residuals between data and total model (in blue). The vertical lines mark all features modeled with pPXF procedure, while the gray bands show the spectral regions we masked during the fitting procedure, corresponding to the gap between the detectors and covering a few emission lines not relevant to this work.

<figure><img src="image_24.png" /><figcaption>Fig. A.1. Examples of continuum and emission lines fitting in a spaxel from the southern nucleus (top panel) and one from the region between nuclei (bottom panel), with three Gaussian components for each line. (See Fig. 3 for details.)</figcaption></figure>

## Appendix B: Moment maps of broad and narrow components

In Figs. B.1 and B.2 we plot the moment maps of [Fe II]λ1.644μm and Paα emission lines. In the upper (lower) panels are shown the maps of the narrow (broad) components.
<figure><img src="image_25.png" /><figcaption>[Fell]λ1.644µm</figcaption></figure>

Fig. B.1. Moment maps of [Fe II]λ1.644µm narrow and broad components. (See Fig. 12 for details.)

<figure><img src="image_26.png" /><figcaption>Paα</figcaption></figure>

Fig. B.2. Moment maps of Paα narrow and broad components. (See Fig. 12 for details.)