© ESO 2022

1. Introduction

Local ultra-luminous infrared galaxies (ULIRGs; with rest-frame [8 − 1000 μm] luminosity, L_IR, in excess of 10¹² L_⊙) are an important class of objects for understanding the formation and evolution of massive galaxies. A classic evolutionary scenario (Sanders et al. 1988; Springel et al. 2005) suggests that ULIRGs evolve into elliptical galaxies through a merger-induced dissipative collapse. In this scenario, the gas of colliding galaxies loses angular momentum and energy, falling into the coalescing centre of the system. Here it serves as fuel for the starburst (SB) and the growth of a supermassive black hole in a dust-enshrouded environment. Then, the system evolves into an optically bright quasar once it either consumes or removes shells of gas and dust through powerful winds. Finally, the merger remnant becomes an elliptical galaxy.

Recent theoretical works have pointed out that dissipative mergers can also lead to the formation of new disk galaxies. Gas that is not efficiently forced to collapse and form new stars, nor expelled by SB and active galactic nucleus (AGN) winds, can be preserved in a disk or re-form a new disk plane and start regrowing a stellar disk (Robertson et al. 2006; Robertson & Bullock 2008; Bullock et al. 2009; Governato et al. 2009; Hopkins et al. 2009, 2013). Hydrodynamical simulations show that cold flows from filamentary structures also play a major role in the buildup of disks in galaxies (Kereŝ et al. 2005; Dekel et al. 2009; Governato et al. 2009). The interaction between inflows and outflows, the amount of gas, and the mass ratio of the merging galaxies and their orbital parameters (e.g. Hopkins et al. 2009) all affect the probability of preserving (or reforming) a disk. From an observational point of view, ordered disk-like kinematics are generally observed in merger systems, both at low z (Bellocchi et al. 2013; Medling et al. 2014; Ueda et al. 2014; Barrera-Ballesteros et al. 2015; Perna et al. 2019) and at high z (up to z ≳ 4; e.g. Hammer et al. 2009; Alaghband-Zadeh et al. 2012; Harrison et al. 2012; Perna et al. 2018; Leung et al. 2019; Tadaki et al. 2020; Cochrane et al. 2021).

Recently, we started a project aimed at studying, at sub-kiloparsec scales, the two-dimensional multi-phase outflow structure in a representative sample of 25 local ULIRGs by comparing the capabilities offered by the Atacama Large Millimeter/submillimeter Array (ALMA) interferometer and the Multi Unit Spectroscopic Explorer (MUSE) at the Very Large Telescope (VLT). The project, Physics of ULIRGs with MUSE and ALMA (PUMA), is described in the first paper of this series, Perna et al. (2021; Paper I hereinafter). First MUSE data results are also presented in Paper I: we derived stellar kinematics for all the PUMA systems and found that post-coalescence systems are more likely associated with disk-like motions, while interacting (binary) systems are dominated by non-ordered and streaming motions. We also investigated the presence of nuclear outflows associated with the individual nuclei and found ionised and neutral outflows in almost all individual nuclei. A more comprehensive study of physical and kinematic properties of the interstellar medium (ISM) of the archetypical ULIRG Arp 220 was presented in Perna et al. (2020, as part of the PUMA project). In Pereira-Santaella et al. (2021; Paper II hereinafter) we instead analysed the ∼220 GHz and CO(2−1) ALMA observations to constrain the hidden energy source of ULIRGs, providing evidence for the ubiquitous presence of obscured AGN that could dominate their IR emission.

The PUMA sample also allows us to investigate the presence of rotational motions in connection with inflows and outflows in dissipative mergers and, therefore, to test the predictions of hydrodynamical simulations. Hence, the present paper is aimed at investigating the prevalence of gas rotational motions in the inner regions of PUMA systems, as well as their (mis)alignments with the stellar component. In this work, we also characterise the kinematic properties of the associated disk structures in terms of inclination, rotational velocity, velocity dispersion, and dynamical mass. Finally, we compare PUMA properties with those of other local (U)LIRGs and high-z populations of normal main sequence (MS) and SB galaxies, studying the variation in the gas velocity dispersion as a function of the star formation rate (SFR) and the ‘starburstiness’ of the system, defined as the ratio between the specific SFR (sSFR; i.e. SFR/M_*) of a galaxy and the sSFR of an MS galaxy with the same z and M_* (δMS = sSFR/sSFR|_MS; e.g. Elbaz et al. 2011).

This paper is organised as follows. In Sect. 2 we briefly summarise the PUMA sample selection and present the data analysis of spectroscopic (MUSE) and photometric (Hubble Space Telescope, HST) data. In Sects. 3.1–3.4 we report the main results obtained from the spectroscopic analysis, in terms of incidence of disk-like motions in the gas and stellar components, and also compare the gas and stellar motions along their kinematic major axes. For those systems with disk-like gas motions, we present in Sect. 3.5 3D-Barolo modelling and infer a kinematic classification in terms of the ratio between rotational velocity and intrinsic velocity dispersion. Section 3.6 presents the study of correlations between the intrinsic velocity dispersion and SFRs and the starburstiness in an extended sample of MS and SB galaxies in the redshift range z ∼ 0.03 − 2.6. Finally, Sect. 4 summarises our conclusions. Throughout this paper, we adopt the cosmological parameters H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7.

2. Methods

2.1. Sample

The PUMA sample is a volume-limited (z < 0.165; d < 800 Mpc) representative sample of 25 local ULIRGs. These targets were selected among the IRAS 1 Jy Survey (Kim & Sanders 1998), the IRAS Revised Bright Galaxy sample (Sanders et al. 2003), and the Duc et al. (1997) catalogue, isolating the sources visible by ALMA and MUSE and uniformly covering the ULIRG luminosity range. The sample was also selected to include an equal number of systems with AGN and SB nuclear activity (based on mid-IR spectroscopy) in the pre- and post-coalescence phases of major mergers (with projected nuclear distances lower than 10 kpc; see Paper I for further details).

So far, we have obtained MUSE observations for 85% of the systems in the sample (21 systems with 31 individual nuclei; see Tables 1 and 2 in Paper I), and ALMA CO(2−1) and ∼220 GHz continuum observations for the entire sample (22 systems, with 32 individual nuclei, have been already presented in Paper II). In this work, we focus on the kinematic properties of the ionised gas of the 20 ULIRGs (with 29 individual nuclei) observed with MUSE, therefore excluding the Arp 220 system whose properties have been extensively described in Perna et al. (2020, but see also Appendix B for a brief description of its kinematics). Information about the MUSE data used in this work are collected in Paper I, Table 2. At the mean distance of our targets (∼400 Mpc), the MUSE spaxel scale, resolution, and field of view (FoV) correspond to ∼0.34 kpc (0.2″), ∼1 kpc (0.6″), and ∼100 × 100 kpc² (60″ × 60″), respectively.

2.2. Data analysis

In this section we describe the MUSE spectroscopic and HST imaging data analysis we followed to characterise the kinematics and dynamics in our PUMA targets.

2.2.1. MUSE spectroscopic analysis

The MUSE data reduction and analysis was executed by following an approach similar to that described in Paper I. We briefly summarise it in the following. The data reduction and exposure combination were carried out by using the ESO pipeline (muse – 2.6.2). The astrometric registration was performed using the Gaia Data Release 2 (DR2) catalogue (Gaia Collaboration 2018) for all but 5 systems, for which we used registered HST optical images as reference (because of the absence of Gaia stars in the MUSE FoV; see Paper I, Sect. 3.1 for more details).

We first fitted and subtracted the stellar continuum from each spaxel. To do so, we initially performed a Voronoi binning (Cappellari & Copin 2003) on the cube to achieve a minimum signal-to-noise ratio (S/N > 16) per bin on the continuum. We then fitted the stellar continuum in each bin through the pPXF code (penalized PiXel-Fitting; Cappellari & Emsellem 2004; Cappellari 2017), using the Indo-U.S. Coudé Feed Spectral Library (Valdes et al. 2004) as stellar spectral templates to model the stellar continuum emission and absorption line systems. We then subtracted the stellar continuum from the total spectra in each spaxel, scaling the fit from bin to each spaxel according to the observed continuum flux (see Paper I, Sect. 5.1).

A slightly different approach was instead used for the two Seyfert 1 in our sample, IZw1 and I01572: in addition to the stellar spectral templates, we made use of an AGN template constructed on the basis of the observed nuclear spectrum, modelled with a combination of a power-law continuum, forbidden, and permitted emission lines as described in Paper I, Appendix B. Because of the point-spread-function blending, this AGN component accounts for a significant fraction of the total emission in the innermost nuclear regions, and rapidly reduces going to radii r ≳ 1″. This step allows us to better reconstruct the stellar velocity field in the nuclear regions with respect to our previous analysis results (see Fig. 5 in Paper I, and Fig. C.1, top-right in this work).

Before proceeding with the fit of the emission lines, we derived a second Voronoi tessellation to achieve a minimum S/N = 8 of the Hα line for each bin. This feature has been preferred to the [OIII]λ5007 line, generally used to trace ionised outflows, as this line is highly absorbed in ULIRG systems due to their large dust content. The use of Hα as a reference for the tessellation allows us to better preserve the important spatial information (both for kinematics and emission line structures).

At this point we fitted the most prominent gas emission lines from the continuum-subtracted cube, by using the Levenberg-Marquardt least-squares fitting code CAP-MPFIT (Cappellari 2017). In particular, we modelled the Hβ and Hα lines, the [O III]λλ4959,5007, [N II]λλ6548,83, [S II]λλ6716,31, and [O I]λλ6300,64 doublets with a simultaneous fitting procedure. To account for broad and asymmetric line profiles, already observed in the nuclear regions of almost all PUMA targets (Paper I), we performed each spectral fit five times at maximum, with one to five kinematic components, that is, Gaussian sets, each centred at a given velocity and with a given full width at half maximum (FWHM). The final number of kinematic components used to model the spectra was derived on the basis of the Bayesian information criterion (BIC; Schwarz 1978). A detailed description of the Gaussian fit routine can be found in Paper I and Perna et al. (2020).

In Fig. 1 we show two examples of our continuum-subtracted, high S/N spectra, extracted from two different regions of the target 13120−5453 (I13120 hereinafter). The best-fit models show the presence of broad and asymmetric line profiles in the two spectra, and a significant diversity in the relative contribution of narrow and broad components: the emission lines in the top a panels are dominated by the contribution of extremely broad Gaussian components, while the ones in the bottom b panels have a well defined narrow core, especially in the Balmer transitions. These two examples show that a multi-component simultaneous fit of all prominent optical emission lines is required to properly separate ordered and perturbed motions in our PUMA systems.

Fig. 1. Examples of our multi-component Gaussian fit decomposition for two continuum-subtracted spectra of I13120, extracted from single spaxels at ∼1″ (top panels) and 2″ (bottom) south-east of the nucleus. The red curves show the best-fit solutions in the Hβ–[OIII] (left) and Hα–[NII] (right) regions; all emission lines are fitted simultaneously. The vertical solid lines mark the rest-frame wavelength of the emission lines, derived from the stellar velocities of the nuclear spectrum; the vertical dot-dashed blue lines mark instead the local stellar systemic (i.e. at the position of the spaxel from which the spectrum is extracted). These examples demonstrate the diversity of emission-line profiles observed in the FoV of a single target: the spectrum in panels a is dominated by broad and blueshifted Gaussian components, while the one in panels b is dominated by bright narrow Gaussians close to the systemic velocity (for all but the [OIII] lines).

2.2.2. Complementary photometric analysis

We made use of the Photutils (Bradley et al. 2016) Isophote package of Astropy (Astropy Collaboration 2018) to perform a basic photometric analysis of ancillary HST near-IR images, available for most of our targets. This analysis provides important morphological parameters to be compared with those inferred from the kinematic analysis described in the next sections.

We fitted a series of isophotal ellipses to each galaxy: Isophote was instructed to hold the centre position constant, whereas the ellipticity (ϵ) and position angle (PA) of the ellipses interpolating the galaxy isophotes were allowed to vary (e.g. Costantin et al. 2017, 2018). Isophote provides the azimuthally averaged surface brightness profile as well as the variation in ϵ and PA as a function of the semi-major axis length. In Table 1 we report the inferred inclinations, derived from ϵ following Willick et al. (1997), and PA at the galactocentric distance of two times the effective radius, R_e, defined as the radius that contains half of the galaxy light (and computed adopting the curve-of-growth method; see e.g. Crespo et al. 2021). Morphological parameters are reported for a small fraction of sources (12 of 29 individual galaxies), as the intrinsic values of PUMA galaxies are usually distorted by merger interactions, and the presence of companion systems. These distortions are also sometimes responsible of significant variations in the morphological parameters at large radii (r > R_e); we therefore report in the table whether the morphological inclination and PA do show significant variations at large distances.

Our rough estimates for R_e are in general agreement with those obtained in previous works, by a factor of 2 (Veilleux et al. 2006; Haan et al. 2011; Kim et al. 2013), for all but the two Seyfert 1 in our sample, IZw1 and I01572. For these two sources we therefore considered the distance obtained by Veilleux et al. (2006), who performed a more rigorous multi-component two-dimensional image decomposition to separate the host galaxy from its bright active nucleus.

In Sect. 3.4 we compare the morphological PAs we derived with Isophote with the kinematic ones derived from MUSE spectroscopic analysis. The inclination measurements will be instead used to model the gas kinematics with 3D-Barolo (Di Teodoro & Fraternali 2015) in Sect. 3.5.

2.3. Emission line tracers and velocity parameters

Throughout this paper, we differentiate between the individual kinematic component parameters, FWHM_j and ΔV_j, with j from 1 to 5 at maximum and the non-parametric velocities v10, v50, v90, and W80 (e.g. Liu et al. 2013). The former identify the width of a specific kinematic component, j, and its velocity shift with respect to the systemic, defined as the stellar velocity in the nuclear position (see Sect. 5.2 in Paper I); FWHM_j and ΔV_j are common to all emission lines fitted simultaneously. Instead, the non-parametric velocities are defined as follows: v10, v50 and v90 are the 10th, 50th, and 90th-percentile velocities, respectively, calculated on the (multi-component) fitted line profile with respect to the systemic. Therefore, they correspond to the velocities at which 10, 50, and 90% of the line flux is accumulated. W80 is defined as v90 − v10.

Observations of relatively large samples of AGN and star-forming galaxies (SFGs) have indicated that the Hα emission is not necessarily dominated by outflows as is the case for [OIII] emission lines (e.g. Bae & Woo 2014; Cicone et al. 2016), as the Balmer line has a significant contribution from star-forming regions. Indeed, the spectral analysis of the line emission coming from the nuclear regions of our PUMA systems already revealed that v10 and W80 of Hα are, on average, 20% smaller than those of the [OIII] (Paper I). Similarly, [NII] lines have a larger velocity width than Hα, consistent with those of the [OIII] (⟨W80 [O III]/W80 [NII] ⟩ ∼ ⟨v10 [O III]/v10 [NII] ⟩ = 1.05). This empirical evidence can be explained by taking into account the fact that [NII] is brighter than Hα both in AGN ionisation cones, often affected by gas flows (e.g. Fischer et al. 2013), and in shocks (e.g. Allen et al. 2008; Perna et al. 2020). Therefore, the nitrogen line may be better analogous to the [OIII] emission, as preferentially traces outflows and more unsettled material compared to Hα (see also e.g. Harrison et al. 2016). Because of the high extinction in PUMA systems, [NII] has to be preferred to the (fainter) [OIII] as outflow tracer (see also Perna et al. 2020).

For each PUMA target, we produced emission line maps for the flux, v50, and W80 non-parametric velocities, obtained considering the total modelled line profiles made up of the sum of the fitted Gaussians. The Hα and [NII] maps of I13120 are presented in Fig. 2. These maps mark the dissimilarities mentioned above: the Hα velocity field appears more regular than the [NII] one; indeed, the Hα map shows slightly smaller velocity widths. The Hα and, in particular, the [NII] velocity maps show a biconical outflow structure, the approaching part to the south and the receding part to the north. Precisely, the v50 map shows high-v [NII] gas blueshifted to the south, within a conical region with a large opening angle, and high-v redshifted [NII] extending to the north. The outflow biconical structure is also associated with high W80 (up to ∼1000 km s⁻¹). We therefore conclude that, with respect to [NII] lines, Hα is less affected by outflows and highly perturbed kinematics.

Fig. 2. I13120 maps. Top: total Hα integrated flux (left), Hα centroid (v50, centre), and line width (W80, right) obtained from the multi-component Gaussian fit. Bottom: similar panels for [NII]; in the W80 map, the red lines indicate the wide biconical outflow along the north-south direction. Masked regions mark the spaxels contaminated by the presence of background and foreground sources and excluded as disturbing the data analysis. The first solid contour is 3σ, and the jump is 0.5 dex. The cross marks the nucleus. North is up, and west is to the right.

In this paper we focus on the disentangling of ionised gas rotation dynamics in the PUMA sample; we therefore present all Hα maps in Appendix A, and leave the [NII] maps to a following investigation.

3. Results and discussion

3.1. Ionised gas kinematic decomposition

The velocity field of I13120 (Fig. 2) shows a regular velocity gradient along the east-west direction in addition to the typical features observed in biconical outflows, such as blue- and redshifted emitting gas with increased line widths in regions preferentially located along the perpendicular direction. In order to understand if this system is rotationally supported, we take advantage of the kinematic decomposition described in Sect. 2.2. Figure 3, left, shows the distribution of ΔV_j and FWHM_j for each Gaussian component j used to model the emission line profiles in I13120. The figure shows a clear trend: highest FWHM_j (≳700 km s⁻¹) are associated with significant blueshifts (ΔV_j ≲ −400 km s⁻¹), while Gaussian components with smaller FWHM_j have |ΔV_j|≲200 km s⁻¹. The highest velocity shifts and line widths are associated with the outflow (see also e.g. Woo et al. 2016 for similar diagrams obtained from SDSS integrated spectra of nearby AGN); instead, the smallest velocities are associated with less perturbed kinematics. To better investigate the presence of rotationally supported motions, we select in the ΔV_j − FWHM_j plane the components with |ΔV_j|< 250 km s⁻¹ and FWHM_j < 400 km s⁻¹ (see e.g. Mingozzi et al. 2019 for a similar approach), and construct a new data cube for the Hα emission, labelled narrow Hα data cube. The line width threshold FWHM_j = 400 km s⁻¹ was chosen taking into account the fact that in our targets the stellar component, more sensitive to gravitational motions, has a velocity dispersion significantly smaller (with σ_* up to 200 km s⁻¹ only in the innermost nuclear regions, due to beam smearing effects). The flux distribution as well as the velocity field and velocity dispersion of narrow Hα are reported in Fig. 3. As expected, the narrow Hα maps display more regular velocity patterns (and significantly lower line widths) with respect to those obtained from the total Hα (and [NII]) lines in Fig. 2.

Fig. 3. I13120 velocity diagram and maps. Left panel: I13120 velocity shift ΔVj − FWHMj diagram (coloured by density in log-scale) for the individual Gaussian components used to model the emission line profiles in the data cube. The dashed red lines isolate the Gaussian components used to reconstruct the narrow Hα data cube, with |ΔVj|< 250 km s−1 and FWHMj < 400 km s−1. Right panels: narrow Hα flux distribution, velocity (v50), and line width (W80) maps, reported in the second, third, and fourth panels, respectively (see Fig. 2 for details).

In the next section we investigate the presence of rotation in our PUMA systems, taking advantage of this kinematic decomposition between more and less perturbed Gaussian components in the ΔV_j − FWHM_j plane, and extracting for each target the position-velocity (PV) diagrams along the kinematic major axis of narrow Hα data cubes.

3.2. Kinematics along the major axis

Figure 4 shows the PV plots along the kinematic major axis position angle (PA^kin) of I13120, for both the total Hα emission (first panel) and the narrow Hα (second panel). A clear velocity gradient from ∼200 km s⁻¹ to ∼ − 200 km s⁻¹ is observed in both panels; therefore, the exclusion of very high velocity components does not introduce or alter significantly this gradient. In the same figure, we report the extracted line velocity centroids and velocity dispersion of the narrow Hα, as well as the stellar velocities (obtained from pPXF analysis; see Paper I) along the stellar component major axis. At this stage, no correction for the beam smearing is performed in the reported velocity dispersion. Both the narrow Hα and stellar components exhibit (i) a well defined velocity gradient along their major axes, and (ii) a peak in the velocity dispersion diagram at the position of the nucleus. These two conditions provide initial evidence for a rotation-dominated system (e.g. Förster Schreiber et al. 2018).

Fig. 4. I13120 PV diagrams along the galaxy major axis. Left two panels: PV maps of the total Hα and the narrow Hα emission. Right panels: extracted line velocity centroids and line width of the narrow Hα, as well as the stellar velocities along the same axis. 1″ corresponds to ∼0.6 kpc at the distance of I13120, as labelled in the third panel.

In Fig. B.1 we show the comparison between the narrow Hα and stellar velocities along their kinematic major axis, for the 18 targets for which we can observe a clear velocity gradient (and measure a PA^kin) together with a peak in the velocity dispersion diagram at the position of the nucleus for at least one component (i.e. gas or stars). The PA^kin measurements, obtained with the python PaFit package (Krajnovic et al. 2006), the velocity amplitudes and median velocity dispersion measured along PA^kin are reported in Table 1 (Cols. 5 to 10), for both stellar and gas components.

The simple visual comparison between gas and stellar kinematics along PA^kin allows us to isolate five nuclei that are reasonably associated with more regular, disk-like kinematics for both gas and stellar components, according to the two criteria highlighted before: IZw1, I10190 W¹, I12072 N, I13120, and I20100 S. The relatively small number of systems with such characteristics is due to the fact that PUMA consists of advanced interacting ULIRGs systems² with nuclear projected separations smaller than 10 kpc (i.e. systems classified as IIIb, IV, and V in the Veilleux et al. 2002 scheme). The small number of sources (five) with stellar and gas disk-like kinematics does not allow us to infer any specific conclusion about the conditions possibly related to more regular motions, also because of the different intrinsic properties of these systems: IZw1 has a small companion at ∼18 kpc; I10190 and I20100 have two nuclei separated by ≳5 kpc, and prominent tidal features; I12072 has two nuclei at a projected distance of ∼2.3 kpc; I13120 has a single nucleus and extended tails and loops surrounding the main body of the galaxy (up to ≳20 kpc from the nucleus; see e.g. Fig. 1 in Privon et al. 2016).

The rest of the sample displays a variety of kinematics. In the last column of Table 1 we distinguish among systems with evidence for disk-like kinematics in gas and stellar components. For instance, I07251 W, I14348 NE and I17208 do not have well defined kinematic properties: they might present disk-like motions, but with kinematic centres possibly not coincident with the nuclear position (see Fig. B.1). We note, however, that these small offsets might also be due to different amount of dust or the presence of tidal streams along the major axis PA (see e.g. Hα flux distribution in Fig. A.14). A better investigation of these offsets is reported in Appendix C. I20087 has peculiar outflow features, reasonably responsible of the observed velocity gradient (with a maximum variation δv_gas ∼ 353 km s⁻¹, a factor of 2.6 higher than δv_*). Finally, there are seven systems with regular stellar kinematics but highly perturbed gas motions, with σ_gas generally higher than 100 km s⁻¹ across the major axis, and without clear trends in the LOS velocity: I00188, I01572, I05189, I14348 SW, I14378, and I19542 (with the possible addition of I16090, affected by poor data quality). A more detailed description of individual targets is reported in Appendix B.

Summarising, we found five systems with regular, disk-like kinematics traced by the narrow Hα and stars on scales of ≳6 kpc (in diameter), IZw1, I10190 W, I12072 N, I13120, and I20100 S, with the possible inclusion of three additional targets (I14348 NE, I17208, and I07251). The remaining targets (21/29 individual nuclei) show more complex gas kinematics, dominated by tidal streams (e.g. I09022 in Fig. A.6), loops (e.g. I05189 in Fig. A.4) and outflows (e.g. I13451 in Fig. A.10) that prevent a clear identification of possible features due to more regular, disk-like motions. On the other hand, the incidence of stellar disk-like motions is slightly higher than rotating gas, with 13 (17, including the more uncertain systems in Table 1) systems out of 29. This higher incidence is probably due to the fact that the stellar component is less affected by non-gravitational perturbations (shocks, outflows). Therefore, the incidence of gas disk-like kinematics in our PUMA sample, of 27% (8 out of 29), has to be considered as a lower limit. In support of this perspective, we note that in a merger the gas has shorter dissipative timescales than stars, and thus it should settle back on a rotating disk earlier than stars (e.g. Springel & Hernquist 2005).

Rodriguez-Gomez et al. (2017) analysed the morphology of ∼18 000 central galaxies at z ∼ 0 from the Illustris cosmological hydrodynamic simulation (Sijacki et al. 2015). They found that, for objects with M_* ≲ 10¹¹ M_⊙, mergers do not seem to play any significant role in determining the galaxy morphology: remnants are associated with both spheroidal and disk-dominated galaxies (see also Sparre & Springel 2017 for similar results). An incidence of ∼27 − 50% for rotating disks in our PUMA sample is therefore consistent with these theoretical predictions.

We stress however that the PUMA sample, with its relatively small number of targets and the different intrinsic properties of each ULIRG, limits the statistical meaning of our results. Among the 8/29 systems with gaseous disk-like kinematics, I10190 W, I14348 N and I20100 SE are associated with less advanced stages of the merger (wide binaries with nuclear separations ≳5 kpc in projection), and we cannot exclude a disk destruction in subsequent phases; IZw1 is instead a minor merger. On the other hand, the four remaining targets have a unique kinematic centre and kiloparsec-scale rotation signatures, regardless the presence of double nuclei (in the binary systems I07251 and I12072) or strong streams (in the remnants I13120 and I17208). Among the other systems, we identified five merger remnants with a stellar disk (I00188, I01572, I05189, I14378, and I19542) but highly perturbed gas kinematics that might prevent the detection of a gaseous disk. These nine targets (four with a gaseous disk and five with a stellar disk) represent the strongest evidence for the preserving (or reforming) of a gaseous disk in major merger processes within the PUMA sample.

3.3. Differences between gas and stellar kinematics along major axis PA^kin

The narrow Hα and stellar PV diagrams in Fig. 4 (third and fourth panels) display significant dissimilarities: the maximum gas and stellar velocity variations along PA^kin are δv_gas ∼ 328 km s⁻¹ and δv_* ∼ 214 km s⁻¹, respectively, while the peak velocity dispersion are σ_gas ∼ 140 km s⁻¹ and σ_* ∼ 210 km s⁻¹. Differences between gas and stellar kinematics along PA^kin are usually observed in (U)LIRGs (e.g. Cazzoli et al. 2014; Crespo et al. 2021), and can be interpreted as due to the presence of different dynamical structures or distinct levels of obscuration (see below).

A precise comparison between stellar and gas rotation along PA^kin can be performed only for a couple of systems in our sample. Among the 8 systems isolated in the previous section (IZw1, I07251 W, I10190 W, I12072 N, I13120, I14348 NE, I17208, and I20100 S), only two targets show regular gas PV diagrams, without significant contributions from perturbed components: I10190 W and I13120³. These two systems show similar maximum stellar and gas velocities, with δv_gas ∼ 200 km s⁻¹ and δv_* ∼ 300 km s⁻¹, and similar maximum gas velocity dispersion, of ∼130 km s⁻¹; vice versa, their maximum stellar velocity dispersion are significantly different, with σ_* ∼ 200 km s⁻¹ in I13120 and ∼80 km s⁻¹ in I10190 W. This difference might be due to intrinsic dissimilarities, as for instance the nuclear obscuration. As dust preferentially obscures young stars, which tend to be dynamically cooler than older stellar populations, an higher obscuration in the nuclear regions could translate into a higher velocity dispersion. The continuum colour of I13120, defined as log (f_R/f_B), with f_R and f_B being the flux at ∼9000 Å and ∼4500 Å, respectively, is a factor of ∼3 higher than I10190 W in the nuclear regions (Paper I). The different R_e of the two systems (see Table 1) could play a role as well: at fixed disk mass, a more compact disk has a steeper inner velocity gradient, resulting into a higher velocity dispersion peak.

Instead, the observed differences between stellar and gas velocity amplitudes in I10190 W and I13120 (δv_gas ∼ 1.5δv_*) can be related to the star formation activity during the merger stages. For instance, numerical simulations by Cox et al. (2006) showed that old stars that are present prior to the merger (i.e. the oldest stellar populations) are the slowest rotators in a merger remnant; on the contrary, younger stars forming during the first passage of the galaxies and the final merger event are the fastest rotators (see their Fig. 7). We might speculate that youngest stars do not significantly contribute to the measured stellar velocity dispersion, as they are more embedded in dusty regions. In Catalán-Torrecilla et al. (in prep.), we will present the stellar population synthesis and its spatial distribution, in order to test this scenario.

3.4. Morpho-kinematic PA (mis)alignments

Figure 5, left, shows the comparison between morphological and (gas and stellar) kinematic major axis PAs, for all PUMA systems where it was possible to determine PA^kin. Following Barrera-Ballesteros et al. (2015), we compare our morpho-kinematic (mis)alignments with those of a control sample of 80 non-interacting galaxies from the CALIFA survey, whose spatial sampling (from ∼0.3 to 1.5 kpc) and FoV coverage (sizes from 7 to 40 kpc) are comparable with those of our galaxies. In the figure, we report the 1:1 relation, with a shaded region that includes 90% of the CALIFA non-interacting sources (i.e. with a misalignment smaller than 22°). The relatively small PUMA sample does not allow us to derive strong conclusions about the general behaviour of ULIRGs systems; nevertheless, we note that 57% (four out of seven) of the PUMA systems have misalignments larger than 22°, and 42% (5/12) have . These results are consistent with those reported by Barrera-Ballesteros et al. (2015), who analysed the morpho-kinematic misalignments in a larger sample of ∼80 interacting CALIFA galaxies, considering both stellar and gas kinematic PAs.

Fig. 5. PA diagrams. Left: comparison between morphological and kinematic major axis PAs for the stellar (yellow) and gas (red) components. The dashed blue line indicates the 1:1 relation; the blue-shaded region (±22°) includes 90% of the morpho-kinematic PA misalignments of a sample of 80 non-interacting galaxies from the CALIFA survey (Barrera-Ballesteros et al. 2015). Right: comparison between stellar and gas kinematic major axis PAs. The dashed line indicates the 1:1 relation; the shaded region (±15°) includes 90% of the kinematic PA misalignments of non-interacting CALIFA galaxies (Barrera-Ballesteros et al. 2015).

Figure 5, right, shows instead the comparison between gas and stellar kinematic PAs, for the 13 PUMA systems where it was possible to determine the major axis PAs. Also in this case, we compare our results with those in Barrera-Ballesteros et al. (2015): in the figure, we report the 1:1 relation, with a shaded region that includes 90% of the CALIFA non-interacting sources (i.e. with a misalignment smaller than 15°). About 38% (5 out of 13) of the PUMA kinematic misalignments are larger than 15°, roughly consistent with Barrera-Ballesteros et al. (2015), who find that 20% of the CALIFA interacting sample has kinematic misalignments larger than 15°.

The most deviating points in Fig. 5, right, are associated with I10190 E, I14348 NE, I16090, I20087, I20100 SE (see Table 1). The slightly larger number of PUMA systems with more extreme kinematic misalignments might be due to their more advanced merger stage with respect to CALIFA interacting galaxies, the sample of which contains ≳43% pre-merger systems without any visual feature of interaction and projected distances of up to 160 kpc. In fact, the presence of more close companions, prominent tidal streams, and strong nuclear winds in our PUMA systems might all contribute to the kinematic misalignments (but see also Chen et al. 2016).

These results indicate that interactions and mergers do have an impact on the internal kinematic alignment of galaxies. However, we note that stellar and gas PAs are roughly aligned (Fig. 5, right), while more significant misalignments can be found between the morphological and kinematic PAs (Fig. 5, left), consistent with the Barrera-Ballesteros et al. (2015) results.

3.5. 3D-Barolo and gas kinematics classification

To test whether the systems with more regular gas kinematics are compatible with a rotationally supported system, we modelled the narrow Hα data cubes with 3D-Barolo (Di Teodoro & Fraternali 2015). In particular, we modelled the gas kinematics of the following systems: IZw1, I07251, I10190 W, I12072 N, I13120, I14348 NE, I17208, and I20100 SE. In this section we present the general strategy adopted for I13120; more details about the fitting procedure per individual targets are reported in Appendix C.

The main assumption of the 3D-Barolo model is that all the emitting material of the galaxy is confined to a geometrically thin disk and its kinematics are dominated by pure rotational motion. The possible presence of residual components associated with the outflow might affect the 3D-Barolo modelling, especially in the innermost nuclear regions, where the outflow is stronger. Nevertheless, this model allows us to assess the presence of such disks and to infer a simple kinematic classification through the standard v_rot/σ₀ ratio, where v_rot is the intrinsic maximum rotation velocity (corrected for inclination, v_rot = v_LOS/sin(i)) and σ₀ is the intrinsic velocity dispersion of the rotating disk, related to its thickness. In this work, we define σ₀ as the measured line width in the outer parts of the galaxy, corrected for the instrumental spectral resolution (e.g. Förster Schreiber et al. 2018).

3D-Barolo best-fit results have been obtained following two different approaches. The first one consists of a two-step strategy. First, we tried different azimuthal models spanning a range of disk inclination angles i with respect to the observer (5 to 85° spaced by 5°, with 0° for face-on); during this step, the i parameter is fixed, and the fitting minimisation is performed considering the following free parameters: v_rot, the rotation velocity, σ, the velocity dispersion, and ϕ, the major axis PA. The disk centre is fixed to the position of the nucleus (inferred from registered HST/F160W images; see Paper I). We therefore inferred the disk inclination angle considering the best-fit configuration with the minimal residuals, defined using the Eqs. (2) and (3b) in Di Teodoro & Fraternali (2015). Then, we run 3D-Barolo with a local normalisation, letting it minimise the v_rot, σ, ϕ and i parameters. In this second step, the inclination is left free to vary in a few degrees around the best-fit i defined in the previous step. For the second method, we simply run 3D-Barolo with a local normalisation, but initialising the inclination to the value derived from the isophote modelling of HST data (Sect. 2.2.2), hence assuming that continuum and narrow Hα have the same geometry. As for the PA measurements, 3D-Barolo fit analysis is performed on the innermost nuclear position, excluding the regions with poor S/N (< 3), for which a Voronoi tessellation would be required.

The resulting best-fit plots for I13120 are shown in Fig. 6, while best-fit parameters are reported in Table 2, together with those of the remaining seven targets with evidence of rotation (see Appendix C). We note that the 3D-Barolo best-fit inclination of I13120 obtained with the two methods, i_I = 34° ±3° and i_II = 27° ±3°, are still consistent with the value we derived from the isophote modelling of HST/F160W data, i_morph = 25° ±11° (Sect. 2.2.2). However, their slightly different i values translate in different rotational velocities; we therefore decided to report in the table the best-fit results obtained with both methods. Similarly, for each target in the table we indicate if both methods provide totally consistent results (I10190 W and I20100 SE) or not (I13120 and IZw1), or alternatively, if the results are obtained from the first (I07251, I12072 and I14348 NE, with no isophote analysis) or the second method only (I17208, with unconstrained i when fitted with the first approach).

Fig. 6. I13120 narrow Hα 3D-Barolo disk kinematic best fit of the moment 0, 1, and 2 (first to third rows) and PV diagrams along the minor and major disk kinematic axes (bottom). The black curves in the velocity dispersion map obtained from the data (third row, left panel) identify the region from which σ0 is extracted. In the PV diagram along the major axis (bottom right) and minor axis (left), data are indicated with a grey-scale map and blue contours, while best-fit models are shown with red contours.

The σ₀ values reported in Table 2 are estimated from the narrow Hα velocity dispersion map, as the median value in a radial elliptical annulus that takes care of the disk inclination (as shown in the velocity dispersion panels; see e.g. Fig. 6). This was preferred to the beam-smearing-corrected value that could be inferred from 3D-Barolo, because of the significant fit residuals in the velocity maps, and for consistency with previous works in the literature (see next sections). The use of an annulus region allows us to mitigate the beam-smearing effects or residual outflow contributions, which are higher in the centre than the outside (see e.g. PV diagrams in Fig. 6), and – more generally – remove different contributions due to tidal streams and companion systems, which would artificially increase the velocity dispersion by ∼20% on average. We used these σ₀ values to measure the ratio v_rot/σ₀ for all the galaxies with indication of gas rotation (Col. 7 in Table 2).

It is important to note that our best-fit models have important limitations and systematic uncertainties, and the small formal errors on a parameter do not necessarily imply a good fit (see Neeleman et al. 2021). Both the very simplified disk models and a (possible imprecise) separation between narrow and perturbed Hα components (Sect. 2.2) might be responsible of the significant residuals we observe in the velocity and velocity dispersion maps (e.g. Fig. 6). Nevertheless, our 3D kinematical analysis shows that, on average, this small sub-sample of PUMA systems have a ratio of rotational velocity to velocity dispersion of v_rot/σ₀ ∼ 1 − 8. Although slightly lower than that of spiral galaxies in the local Universe (v_rot/σ₀ ∼ 10), our values are still comparable with Hα measurements of other low-z (U)LIRGs in the literature (e.g. Bellocchi et al. 2013; Crespo et al. 2021) and systems at z ∼ 0.5 − 1 (Rizzo et al. 2021 and references therein). Therefore, 3D-Barolo results provide further indication of rotationally supported gas motions in these targets.

3.6. Gas velocity dispersion in (U)LIRGs and high-z populations: Dependence on starburstiness

Many theoretical and observational studies suggest that gas in high-z galaxies has larger random motions compared to nearby galaxies: in particular, the ionised gas velocity dispersion goes from ∼20 km s⁻¹ in nearby spirals to ∼45 km s⁻¹ in massive MS star-forming disk galaxies at z ∼ 2, although with a significant scattering of values (e.g. Übler et al. 2019).

Figure 7, top, shows the narrow Hα velocity dispersion of our rotationally supported PUMA systems (red circles) as a function of the redshift, together with the Übler et al. (2019) evolutionary trend of SFGs. This trend mostly traces the velocity dispersion evolution of normal MS galaxies (e.g. Übler et al. 2019; Förster Schreiber et al. 2018); we therefore labelled it as σ_0, MS hereinafter. The σ₀ of our PUMA systems, in the range 30 − 85 km s⁻¹, are more compatible with – or possibly higher than – those of high-z galaxies rather than nearby spirals. This can be explained taking into account the following arguments. On the one hand, the velocity dispersion increases as natural consequence of the availability of huge gas reservoirs and intense star formation that is taking place in ULIRGs and high-z galaxies (e.g. Lehnert et al. 2009; Arribas et al. 2014; Johnson et al. 2018). On the other hand, both the gravitational instabilities due to the galaxy interactions and the presence of non-circular motions in our PUMA targets can contribute to the σ₀ enhancement with respect to isolated nearby galaxies.

Fig. 7. Gas velocity dispersion σ0 evolutionary trend and correlations. Top: velocity dispersion, σ0, as a function of z for the PUMA sub-sample (red dots) and other individual ionised gas measurements from the literature, distinguishing between the different samples presented in Appendix D, as labelled. The Übler et al. (2019) evolutionary trend of MS galaxies is shown with a solid curve (shaded area: 1σ scatter around the average trend). Middle: σ0 as a function of the SFR for all targets already reported in the first panel and the KMOS3D galaxies (grey points; Übler et al. 2019). K18 models and a linear fit (with scatter) are also reported, as labelled. Bottom: σ0 normalised to the evolutionary trend of MS galaxies as a function of δMS for all targets already reported above. A linear fit is also reported.

To better investigate the origin of the differences between PUMA systems and normal MS galaxies at different redshifts, in Fig. 7 we report additional individual ionised gas measurements of SB disk galaxies from the literature (see Appendix D for details). Many of them show a significant deviation from the σ_0, MS evolutionary trend, similar to our PUMA systems.

All these measurements from the literature have been obtained from integral field spectrograph data; therefore, they are not strongly affected by beam-smearing effects and other systematics that tend to overestimate the intrinsic dispersion (see discussion in Übler et al. 2019). We also stress here that all selected individual sources are disk galaxies. This ensures relatively small contribution of outflows in the velocity dispersion measurements, the incidence of which increases with SFR and AGN activity (e.g. Cicone et al. 2016; Villar Martín et al. 2020). A few additional caveats should be kept in mind regarding this compilation of sources: all of them are presented as SB galaxies in the original papers, but we note that (i) there is no rigorous definition of a SB galaxy, but several different criteria are often used, and (ii) especially at high z, stellar mass and SFR measurements can be highly uncertain, depending on the availability of multi-wavelength information. This aspect is further discussed below.

3.6.1. The σ₀–SFR correlation

Most of the galaxies presented in Fig. 7, top, significantly deviate from the σ_0, MS evolutionary trend. In order to understand if this deviation is due to the extreme SFR in these systems, we show in Fig. 7, middle, the velocity dispersion σ₀ as a function of the SFR, for all the targets already mentioned, in addition to the sample of normal MS galaxies used to derive the σ_0, MS trend (Übler et al. 2019). All SFR measurements of SB reported in the figure are obtained from IR luminosities (as reported in the original papers, or using the Kennicutt 1998 relation, assuming the Chabrier initial mass function). For the PUMA targets, we considered the nuclear IR luminosities reported in Paper II, Table 7; for the binary systems, the fraction of the IR luminosity assigned to each nucleus is based on their relative ALMA continuum fluxes. We observe a relatively poor correlation (Spearman rank correlation coefficient 0.4), reasonably due to a bias selection: a clearer correlation is in fact observed when combining samples of galaxies that cover a larger dynamical range (i.e. also including local MS galaxies; e.g. Yu et al. 2019; Varidel et al. 2020). By performing a linear regression fit we derive

(1)

which is compatible (within 1σ) with the results from Arribas et al. (2014), obtained from a sample of (U)LIRGs observed with the optical spectrograph VIMOS. We compare this fit with the model predictions by Krumholz et al. (2018), for both high-z galaxies and local ULIRGs (dash-dotted lines in Fig. 7, middle). The slight differences between the two theoretical curves in the figure are due to the distinct ISM physical conditions of these two classes of sources (see Table 3 in Krumholz et al. 2018). These models explain the observed slow increase in σ₀ as a function of SFR in the range [10⁻³, 10] M_⊙ yr⁻¹, followed by a steeper increase up to 10s km s⁻¹ considering two different regimes. The σ₀ floor is due to stellar feedback processes, while at higher SFR the velocity dispersion is regulated by gravitational turbulence. We note that the comparison between these theoretical curves and our collected data is strongly limited by caveats. To begin, all SB galaxies in our plot are disk galaxies: this could result in the exclusion of all targets with more extreme (i.e. higher) σ₀. The next problem is that the Bellocchi et al. (2013) and PUMA velocity dispersion measurements in this plot are derived from the narrow Hα (i.e. after removing the more extreme kinematic components). As shown in Figs. A.1–A.19 and B, significantly higher velocity dispersion (up to several 100s km s⁻¹) is measured in the total Hα line profiles of our PUMA galaxies, reasonably associated with extended outflows and streaming motions. Finally, the ISM of (U)LIRGs and high-z SB galaxies might not be in vertical pressure nor energy balance, as assumed in Krumholz et al. (2018) models: instead, their velocity dispersion might be strongly affected by non-circular motions. Specifically, these σ₀ measurements might not be dominated by the turbulent component, which is relevant for a comparison with the SFR (e.g. Bacchini et al. 2020).

3.6.2. The σ₀/σ_0, MS − δMS correlation

In the σ₀–SFR plane, local (U)LIRGs and high-z galaxies tend to occupy the same region, regardless of their intrinsic difference in terms of morphology, gas fraction and starburstiness. In order to distinguish between normal MS and SB galaxies, we show in Fig. 7, bottom, the velocity dispersion σ₀ normalised to σ_0, MS (solid line in the top panel; Übler et al. 2019) as a function of the starburstiness δMS = sSFR/sSFR|_MS for all targets already reported in the previous panels. The sSFR|_MS is derived from the Speagle et al. (2014) relation, starting from the available stellar mass measurements from the literature for Johnson et al. (2016), Förster Schreiber et al. (2018), Molina et al. (2020), Cochrane et al. (2021), and KMOS^3D individual targets; stellar masses of LIRGs and ULIRGs from this study, from Bellocchi et al. (2013), Pereira-Santaella et al. (2019) and Crespo et al. (2021), are instead derived from the dynamical mass estimates, assuming M_* = (1 − f_gas)×(1 − f_DM)×M_dyn, where f_gas and f_DM are the gas and dark matter fractions⁴, respectively, and M_dyn is the dynamical mass within 2R_e (see next section). For the gas fraction, we considered a conservative f_gas = 0.1 (Isbell et al. 2018; higher f_gas values would further increase their δMS). This gas fraction is consistent with the estimate we obtain considering the molecular gas mass inferred by ALMA data (Paper II; Lamperti et al., in prep.) and the M_* measurements available for a few PUMA targets (see Table 2), . For the dark matter fraction (within 2R_e) we assumed f_DM = 0.26, defined as 1 − M_bar/M_dyn, with M_bar = M_* + M_gas = M_*/(1 − f_gas). The f_DM estimate was derived for the PUMA and Johnson et al. (2016) SB galaxies for which stellar masses are available, and considering the dynamical masses within 2R_e (see Sect. 3.7). Because of these assumptions, we considered a factor of 3 uncertainties for the stellar mass measurements of (U)LIRGs. These uncertainties, however, play a minor role in the derived δMS: at low z, the MS has a soft slope, and normal and massive MS galaxies have similar SFRs (e.g. at z ∼ 0.1, galaxies of 10¹⁰ and 10¹¹ M_⊙ have SFR|_MS ∼ 1 and ∼3.5 M_⊙ yr⁻¹, respectively); on the other hand, local (U)LIRGs have much higher SFRs, from 10 s to 100 s M_⊙ yr⁻¹, and therefore δMS of the order of 10 − 100. Finally, for the Alaghband-Zadeh et al. (2012) and Harrison et al. (2012) high-z galaxies we assumed that M_* = 10¹¹ M_⊙, following Harrison et al. (2012).

Figure 7, bottom, shows a clear correlation between σ₀/σ_0, MS and δMS (Spearman rank correlation coefficient 0.6), suggesting that SB galaxies tend to have higher velocity dispersion than normal galaxies at given z and stellar mass. Tentative evidence of such correlation was already reported by Wisnioski & Förster (2015), for the KMOS^3D galaxies at z ∼ 2, and by Varidel et al. (2020), for nearby MS galaxies of the SAMI survey, but the lack of dynamical range in terms of starburstiness in these surveys maintained the correlation at low significance level (see also Figs. 15 and 16 in Übler et al. 2019).

By performing a linear regression fit, we derived

(2)

Given the correlation between σ₀ and SFR (e.g. Varidel et al. 2020) and the tight inter-relationship between δMS and f_gas (e.g. Wisnioski & Förster 2015; Tacconi et al. 2020), it is unsurprising that the starburstiness correlates with the excess in the velocity dispersion with respect to MS galaxies. The positive slope of ∼0.27 is inconsistent with the one observed between f_gas and δMS, +0.5 (Tacconi et al. 2020), suggesting that complex interactions between different physical drivers are responsible of the correlation observed in Fig. 7 (bottom). As a final check, we studied the correlation between σ₀ and v_rot, as well as between σ₀/σ_0, MS and v_rot, obtaining Spearman coefficients of ∼ − 0.2 (with p-value ∼ 0.05); this shows that not even v_rot can have a significant role in determining the increase in the velocity dispersion of SB galaxies.

The strong correlation σ₀/σ_0, MS − δMS suggests that the SFR of galaxies above the MS is taking place in an ISM significantly more unsettled than in normal (i.e. MS) galaxies. This is likely due to the presence of interactions and mergers, which enhance SFRs while simultaneously increasing the velocity dispersion of the ISM. The absence of a strong correlation between the SFR and the elevated velocity dispersion in star-forming clumps both in local (U)LIRGs and high-z SFGs (Arribas et al. 2014; Genzel et al. 2011) further suggests that the extreme dispersion cannot simply be related to the strong SFR in these systems.

Our arguments are consistent with the parsec-resolution hydrodynamical simulations of major mergers presented by Renaud et al. (2014): they found that the increase in ISM velocity dispersion precedes the star formation episodes. Therefore, this enhancement is not a consequence of stellar feedback but instead has a gravitational origin. In this scenario, δMS can be interpreted as a tracer of the strength of gravitational torques: stronger gravitational torques during the interactions lead the gas to flow inwards, both increasing the velocity dispersion and the efficiency in converting gas into stars. We finally mention that AGN outflows, which are ubiquitous in these systems (Paper I; Paper II), can also contribute to increasing the velocity dispersion of the disks, as suggested by high-resolution hydrodynamic simulations (e.g. Wagner et al. 2013; Cielo et al. 2018).

A more detailed investigation of the physical meaning of the correlations reported in Fig. 7 goes beyond the purpose of this study; here we just stress that, by selecting a (relatively small) sample of MS and SB disk galaxies in the redshift range 0.03 − 2.6, a more significant correlation is observed between σ₀/σ_0, MS and δMS rather than between σ₀ and SFR. We argue that this result might be even more evident considering the entire population of (U)LIRGs (i.e. without excluding targets with no evidence of a rotating disk; see e.g. Table 1), and measuring the velocity dispersion without excluding possible contribution from outflows and streaming motions (see Figs. A.1–A.19).

3.7. Dynamical masses

In this section we derive the dynamical masses of our PUMA sub-sample, and compare them with those of other (U)LIRGs from the literature. Assuming that the source of the gravitational potential is spherically distributed, we can estimate the dynamical mass within a radius R as:

(3)

where G is the gravitational constant, v_circ is the circular velocity in km s⁻¹, and R is given in kpc.

We used the near-IR continuum 2R_e as the radius to calculate M_dyn, which for an exponential profile contains 85% of the total flux. The effective radius of I10190 W, I13120 and I17208 is computed with the Isophote package of Astropy (Sect. 2.2.2); IZw1 R_e is instead taken from Veilleux et al. (2006), who performed multi-component two-dimensional image decomposition to separate the host galaxy from its bright active nucleus; for the remaining three targets, we adopted the ULIRGs average R_e derived by Bellocchi et al. (2013), as the presence of nearby nuclei (I07251 and I12072) and strong tidal features (I14348 NE) do not allow us to model the continuum with isophotal ellipses.

To infer the circular velocity we consider both the rotation and dispersion motions traced by the narrow Hα. In particular, we included the asymmetric drift term, which represents an extra component due to the dispersion of the gas around the disk of the galaxy,

(4)

The term η is a constant and can vary between approximately 1.5 and 6, depending on the mass distribution and kinematics of the galaxy: higher values indicate higher turbulence in the ISM of a rotating disk (Neeleman et al. 2021). We assumed η = 3, following Dasyra et al. (2006), which is very close to the value expected for an exponential, turbulent pressure-supported disk (η = 3.4), and considered an uncertainty of 1.5 to take into account the large range of possible values. We note that, on average, our dynamical masses would be a factor of 1.15 lower (1.32 higher) assuming η = 1.5 (6).

The assumption of a spherically distributed ISM is at odds with that we used in Sect. 3.5 to measure the rotational velocities in our systems. In order to account for this, following Neeleman et al. (2021), we conservatively increased the dynamical mass uncertainty by 20% towards lower masses. This corresponds to consider that the effective total mass distribution falls somewhere in between a thin disk and a sphere.

The measured dynamical masses of our PUMA systems range from ∼2 to ∼7 × 10¹⁰ M_⊙, consistent with the median value derived by Bellocchi et al. (2013) for ULIRG systems, 4.8 × 10¹⁰ M_⊙, confirming that ULIRGs are intermediate mass systems, as previously suggested (i.e. Colina et al. 2005; Rodríguez Zaurín et al. 2010).

In this final part, we further discuss about the dark matter fraction reported in the previous section, and defined as f_DM = 1 − M_bar/M_dyn. Using a sub-sample of 27 SB disk galaxies with available M_* measurements, and assuming f_gas = 0.1, we obtained a median value f_DM = 0.26. We note, however, that for a few sources (6/27, mostly from the PUMA sample) M_* > M_dyn (Table 2; see also Table 6 in Rodríguez Zaurín et al. 2010). This can suggest either that the systems are not relaxed (due to the interaction and mergers) and that M_dyn is unreliable, or that M_* measurements have high uncertainties. As the M_dyn estimates are in agreement with previous works, and because of the fact that for binary PUMA systems in Table 2 the available M_* measurements are obtained without separating the contribution of the merging galaxies, we favour the second interpretation: that M_* measurements are highly uncertain, though a combination or both may also be possible. These arguments led us to consider significant (factor of 3) uncertainties in the determination of the stellar masses for the entire sample of local (U)LIRGs, required to estimate their δMS. The conclusions reported in the previous section are however not affected by these uncertainties, because of extreme SFR in (U)LIRGs targets.

4. Summary and conclusions

The PUMA project is a survey of 25 nearby ULIRGs observed with MUSE and ALMA. It is a representative sample that covers the entire ULIRG luminosity range and includes a combination of systems with AGN and SB nuclear activity in advanced interacting and merging stages. Paper I presents the first MUSE results on the spatially resolved stellar kinematics and the incidence of ionised outflows in nuclear spectra; Paper II analyses high-resolution (400 pc) ∼220 GHz continuum and CO(2−1) ALMA observations to constrain the hidden energy sources of ULIRGs. In this paper, we have investigated the presence of rotational ionised gas dynamics in PUMA targets to understand if, as predicted by models, rotation disks can be preserved during the merging process (or rapidly regrown after coalescence) and, if so, what their main properties are. Our results are summarised below.

(a) We have presented the spatially resolved Hα flux and kinematic maps for the entire PUMA sample, obtained from multi-component Gaussian fit analysis (Figs. A.1–A.19). Irregular large-scale ionised gas velocity fields associated with tidally induced motions and outflows are found in almost all targets; Hα velocities (v50) of up to ∼ ± 300 km s⁻¹ are detected in the MUSE FoV, while HαW80 line widths range from ∼100 to ∼1500 km s⁻¹. The [NII] (and [OIII]) line transitions are even more affected by perturbed motions, such as tidal streams and outflows.

(b) We have studied the Hα kinematics to infer the presence of rotating disk signatures. A kinematic decomposition was performed by selecting in the ΔV_j − FWHM_j plane all best-fit Gaussian components with relatively small velocities and constructing new narrow Hα data cubes. In these newly generated data cubes, the emission associated with gas components with extreme velocities (likely due to outflows and/or tidally driven flows) is minimised.

(c) By studying the gas kinematics along the major axes of our galaxies in the innermost regions (∼5 − 20 kpc), we found that 27% (8 out of 29) individual nuclei are associated with disk-like motions. This has to be considered as a lower limit as the presence of vigorous winds and gravitational torques, as well as observational limitations (in terms of spatial resolution, spectral resolution, and S/N), limits our capabilities to isolate more regular, disk-like kinematics through a multi-component Gaussian fit decomposition. This is supported by the fact that five merger remnants in our sample present stellar disk motions but highly perturbed gas kinematics. Indeed, the incidence of rotating ionised gas disks is a factor of ≲2 smaller than that of stellar disk-like motions (Paper I). This possibly suggests that (i) we are actually missing a significant fraction of sources with gas rotation because of the above-mentioned limitations, or (ii) the gas component is more affected by winds and gravitational interactions, and the probability of preserving a gas disk is lower than that of a stellar disk. In both instances, our results show that, as predicted by models, rotation disks can be preserved during the merging process and/or rapidly regrown after coalescence.

(d) For the eight galaxies with evidence of disk-like motions, we modelled the narrow Hα data cubes with 3D-Barolo and derived rotational velocities v_rot ∈ [70 − 300] km s⁻¹. By combining them with the measured velocity dispersion, σ₀ (∈[30 − 80] km s⁻¹), we derive v_rot/σ₀ values in the range 1−8, providing a further indication of rotationally supported gas motions in these ULIRGs. We also derived their M_dyn, obtaining values in the range (2 − 7)×10¹⁰ M_⊙, consistent with the M_dyn of other ULIRGs in the literature.

(e) We compared the narrow Hα velocity dispersion, σ₀, of our eight PUMA disk galaxies with those of other SB and normal MS disk galaxies at low and high z. We found that all SB galaxies tend towards higher σ₀ values compared to MS galaxies at the same redshift. Interestingly, when we normalise σ₀ to the value expected for MS galaxies (at the same z), considering the Übler et al. (2019) evolutionary trend, σ_0, MS, we found a significant correlation between σ₀/σ_0, MS and the starburstiness, δMS. In particular, SB galaxies display velocity dispersions that are up to a factor of ∼4 higher than those of normal MS galaxies at the same redshift. The relatively poor correlation between σ₀ and the SFR (Fig. 7, middle) suggests that stellar activity cannot be the main mechanism responsible for the σ₀ enhancement observed in SB galaxies, and other mechanisms, possibly related to interactions and mergers, should be taken into account (see e.g. Renaud et al. 2014).

We note, however, that most of the SB galaxies at z ≳ 0.4 collected from the literature are consistent with δMS = 1 once homogeneous recipes are used to derive the SFR and measurement uncertainties are taken into account. As a result, the correlation reported in the figure is mostly driven by the comparison between z ∼ 0.03 − 0.4 (U)LIRGs and KMOS^3D MS galaxies at z ∼ 0.6 − 2.6. This makes a further investigation of gas dynamical conditions in SB galaxies at z > 0.4 highly desirable. The James Webb Space Telescope NIRSpec integral field spectrograph, with its wide spectral range (from 0.6 to 5.3 μm) and sub-arcsecond resolution, will allow a comprehensive characterisation of the ionised gas dynamical conditions in such systems.

Acknowledgments

We thank the referee for an expert review of our paper. The authors thanks Elena Valenti for her support when preparing the observations, and Giustina Vietri for useful discussion on spectral analysis of type 1 AGN. M.P. is supported by the Programa Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715. M.P., S.A., C.T.C. and L.C. acknowledge support from the Spanish Ministerio de Economía y Competitividad through the grant ESP2017-83197-P, and PID2019-106280GB-I00. M.P.S. and I.L. acknowledge support from the Comunidad de Madrid through the Atracción de Talento Investigador Grant 2018-T1/TIC-11035 and PID2019-105423GA-I00 (MCIU/AEI/FEDER,UE). H.Ü. gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. L.C. acknowledges financial support from Comunidad de Madrid under Atracción de Talento grant 2018-T2/TIC-11612 and the Spanish Ministerio de Ciencia, Innovación y Universidades through grant PGC2018-093499-B-I00. E.B. acknowledges support from Comunidad de Madrid through the Attracción de Talento grant 2017-T1/TIC-5213. S.C. acknowledge financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709). A.C.G. acknowledges support from the Spanish Ministerio de Economía y Competitividad through the grant BES-2016-078214. R.M. acknowledges ERC Advanced Grant 695671 “QUENCH” and support by the Science and Technology Facilities Council (STFC). J.P.L. acknowledges financial support by the Spanish MICINN under grant AYA2017-85170-R.

References

Appendix A: Multi-component Gaussian fit results

Fig. A.1. I00188 maps: Hα integrated flux (left), Hα centroid (v50, centre), and line width (W80, right) obtained from the multi-component Gaussian fit. The first solid contour is 3σ, and the jump is 0.5 dex. The cross marks the nucleus. North is up, and west is to the right.

Fig. A.2. IZw1 maps. See Fig. A.1 for details.

Fig. A.3. I01572 maps. See Fig. A.1 for details.

Fig. A.4. I0589 maps. See Fig. A.1 for details.

Fig. A.5. I07251 maps. See Fig. A.1 for details.

Fig. A.6. I09022 maps. See Fig. A.1 for details.

Fig. A.7. I10190 maps. See Fig. A.1 for details.

Fig. A.8. I11095 maps. See Fig. A.1 for details.

Fig. A.9. I12072 maps. See Fig. A.1 for details.

Fig. A.10. I13451 maps. See Fig. A.1 for details.

Fig. A.11. I14348 maps. See Fig. A.1 for details.

Fig. A.12. I14378 maps. See Fig. A.1 for details.

Fig. A.13. I16090 maps. See Fig. A.1 for details.

Fig. A.14. I17208 maps. See Fig. A.1 for details.

Fig. A.15. I19297 maps. See Fig. A.1 for details.

Fig. A.16. I19542 maps. See Fig. A.1 for details.

Fig. A.17. I20087 maps. See Fig. A.1 for details.

Fig. A.18. I20100 maps. See Fig. A.1 for details.

Fig. A.19. I22491 maps. See Fig. A.1 for details.

Appendix B: Position-velocity diagrams

In Fig. B.1 we show the comparison between the (total and narrow) Hα and stellar velocities along their kinematic major axis, for all targets for which we can observe a clear velocity gradient (and measure a PA^kin) together with a peak in the velocity dispersion diagram at the position of the nucleus for at least one component (i.e. gas or stars). The nuclear positions are inferred from registered HST/F160W images (Paper I). The PV plots show only a small portion of the total extension of the ULIRG systems, in order to exclude the contribution from tidal tails, extended outflows or second nuclei, and better identify rotation-like signatures.

Fig. B.1. Position-velocity diagrams along the galaxy major axis for all PUMA systems that show i) a well-defined velocity gradient along the major axis and ii) a peak in the velocity dispersion diagram close to the position of the nucleus. These two conditions provide initial evidence for rotation-dominated kinematics. All systems whose gas (stellar) kinematics satisfy the two conditions are marked with a blue (red) check mark in the top-right corner of the velocity dispersion panel; on the contrary, the systems that do not satisfy at least one of the two conditions are marked with a cross symbol; more uncertain kinematics are marked with a question mark. The systems for which neither the gas nor the stellar components satisfy the two conditions are not reported in this figure. Line velocity centroids and line widths of narrow Hα, total Hα, and stars are reported for each target, as labelled. See Fig. 4 for further details about the individual panels.

Fig. B.1. continued.

The major axis PA^kin measurements have been obtained with the python PaFit package (Krajnovic et al. 2006), for both stellar and gas kinematics (Cols. 5 and 6 in Table 1). The difficulty in measuring reliable PA^kin and obtaining regular velocity profiles led to the exclusion of the following systems: I07251 E and W, I09022, I11095, I12072 S, I13451 E and W, 19297 S, 20100 NW, and I22491 E and W.

For the sources for which the gas and stellar major axes agree within the errors, the PV diagrams have been extracted considering a PA^kin equal to the weighted mean of PA and PA, with the weighting factors equal to the inverse of the quadratic uncertainty on the PA^kin measurements; for those sources for which no PA^kin measurement can be obtained for the gas, we considered the PA (I01572, I05189, I19297 N, I19542); vice versa, PA has been chosen when the stars do not show a clear velocity gradient (I07251 E, I22491 E). Finally, for those sources with different stellar and gas kinematics (I10190 E, I14348 NE, I16090, I20087 and I20100 SE) we extracted the velocity profiles along different PAs.

Before briefly discussing the PV diagrams of individual targets shown in the figure, we note that a clear but shallow velocity gradient (i.e. with δv ≪ 100 km s⁻¹) could translate in a flat σ profile under the MUSE observing conditions (i.e. with an angular resolution ≈0.6″). This would determine the exclusion of disk-like candidates, according to two selection criteria mentioned at the beginning of this section. However, most of PUMA individual systems show δv ≳ 100 km s⁻¹ (see Table 1); the few systems with lower δv (in particular, the gas component in I00188, I22491 E, I14348 SW and I14378) also present disturbed kinematics, excluding the presence of disks even in these conditions.

B.1. Notes on individual targets in Fig. B.1

I00188 shows a regular gradient in V_* and a peak in σ_* at the position of the nucleus, as expected for rotation-dominated kinematics; on the other hand, the narrow Hα velocity profiles is irregular, while its velocity dispersion is close to ∼110 km/s along the entire extension of the major axis, indicating highly perturbed gas kinematics.

IZw1 The Hα in the central pixel is saturated; this translates in a gap in the PV slices shown in the figure. The presence of strong Sy1 emission also prevents us from correctly inferring proper V_* and σ_* measurements in the vicinity of the nucleus; in particular, the σ_* values might be strongly overestimated.

I01572 As for IZw1, the presence of a strong Sy1 prevents us from properly inferring stellar velocity and velocity dispersion in the vicinity of the nucleus. Moreover, the gas kinematics are strongly disturbed by the nuclear outflow and tidal motions.

I05189 shows a clear gradient in V_* and a peak in σ_* at the position of the nucleus, as expected for rotation-dominated kinematics; on the other hand, the narrow Hα velocity profiles is irregular, while its velocity dispersion is close to ∼110 km/s along the entire extension of the major axis, indicating highly perturbed gas kinematics.

I07251 presents disk-like motions, but with a kinematic centre not coincident with the position of either the two nuclei of this interacting system (see also Fig. C.2, top panels). The PV diagrams presented in Fig. B have been therefore extracted considering the kinematic centre obtained from the 3D-Barolo analysis.

I10190 This system shows two rotating disks associated with the two nuclei, separated by ∼7.2 kpc. However, both the stellar and gas kinematics in the vicinity of the E nucleus are strongly affected by the presence of the W system motions.

I12072 This source presents two nuclei, at a projected distance of 2.3 kpc. Disk-like motions are presumably associated with the N nucleus. Broadly regular PV diagrams are observed for both stellar and gas components.

I13120 This target is extensively presented in the main text of this work.

I14348 is an interacting ULIRG with two nuclei, at a projected distance of ∼5 kpc. The NE nucleus presents broadly regular PV diagrams for both gas and stellar components; the SW nucleus presents broadly regular stellar kinematics, and more irregular gas motions due to the presence of a strong outflow.

I14378 shows a clear gradient in V_* and a peak in σ_* at the position of the nucleus, as expected for rotation-dominated kinematics; on the other hand, the narrow Hα velocity profiles is irregular, and σ_gas ≈ 90 km/s along the entire extension of the major axis indicates highly perturbed gas kinematics.

Arp220 shows a broadly regular velocity gradient in both gas and stellar components (see also Scoville et al. 1997). However, the presence of the two nuclei with distinct rotation features on scales of ∼100 pc (at a distance of ∼370 pc), and a kiloparsec-scale, wide-angle outflow prevents a detailed characterisation of the gas kinematics in this system (see e.g. Fig. 19 in Perna et al. 2020). For completeness, in Fig. B.1 we report the PV plots extracted along a PA^kin ∼ 40°, with respect to the position of the two nuclei.

I16090 shows broadly regular gradient in V_* and V_gas, and a peak in σ_* at the position of the nucleus, as expected for rotation-dominated kinematics; on the other hand, σ_gas ≈ 130 km/s along the entire extension of the major axis indicates highly perturbed gas kinematics.

I17208 shows broadly regular velocity profiles in both stellar and gas components; the velocity dispersion profiles are instead more irregular, with higher values towards the south-east direction, probably due to streaming motions (see also Fig. C.6).

I19297 No evidence of disk-like motions is present in the vicinity of the two nuclei of this system. The very high gas velocity dispersion (> 100 km/s) suggests the presence of strongly disturbed kinematics.

I19542 shows a clear gradient in V_* and a peak in σ_* at the position of the nucleus, as expected for rotation-dominated kinematics. On the other hand, the narrow Hα velocity profiles is irregular, mostly associated with blueshifted emission, and the σ_Hα ≈ 110 km/s along the entire extension of the major axis indicates highly perturbed gas kinematics.

I20087 shows a clear gradient in V_* and a peak in σ_* at the position of the nucleus, as expected for rotation-dominated kinematics. The gas velocities broadly resemble the V_* profile, but reaching maximum velocities at ∼2 kpc from the nuclear position a factor of ∼2.6 higher than V_*. This behaviour might be due to the presence of a bi-conical outflow, and will be better investigated in a forthcoming paper.

20100 SE shows broadly regular velocity profiles in both stellar and gas components. We however observe a significant misalignment between the gas and stellar major axis PAs.

Appendix C: 3D-Barolo analysis

As already reported in Sect. 3.5, we followed two different methods to derive 3D-Barolo best-fit results. With the first method, we first tried to constrain the disk inclination using different azimuthal models spanning the almost entire range of inclinations (i ∈ 5 − 85°), selecting the i value with minimal residual. Finally, we ran 3D-Barolo with local models, fitting the rotation velocity v_rot, the velocity dispersion σ, and the major axis PA ϕ, using the disk inclination derived in the first step as initial guess. With the second method, we directly fit all disk kinematic parameters with a local model, by initialising the inclination to the value derived from the Isophote modelling of HST data (Sect. 2.2.2). The general fitting procedure is here explained in more details for each target in our sub-sample of PUMA systems with evidence of rotation.

IZw1 In this target, a correct ionised gas kinematic decomposition between disk and outflow components in the innermost nuclear regions is challenging. The strong outflow, the BLR and iron emission close to the Hα are responsible of a strong degeneracy in the multi-component Gaussian fit results. As a consequence, the narrow Hα map shown in the top-left panel of Fig. C.1 display a blueshifted kinematic component in the innermost nuclear regions not associated with disk kinematics. Similarly, the more extreme redshifted velocities in the south-east direction at ∼3″ are reasonably due to the same fit degeneracy (see also Fig. A.2). On the contrary, the stellar kinematics (top-right panel in Fig. C.1) are more regular and display a clear rotation pattern.

Fig. C.1. IZw1 velocity maps and the 3D-Barolo disk kinematic best fit. Top panels: Narrow Hα (left) and pPXF stellar (right) velocity maps. The magenta line identifies the major axis PA measurement, computed within the black box region; the black cross identifies the nucleus. Second to fourth row panels: Comparison of 3D-Barolo data and model moment maps, as labelled in the figure. In the intensity and velocity maps, we report the major axis PA (dashed line) and the position of the kinematic centre (black cross); the green curves in the velocity maps represent the zero-velocity axis. In the dispersion map on the left, the black curves show the region from which the median σ0 value has been derived. Bottom panels: PV diagram along the major axis (right) and minor axis (left) for both the data (grey map and blue contours) and the best-fit model (red contours; the yellow dots associated with individual concentric rings are used to model the data).

The results obtained with I and II methods are broadly consistent (see Table 2), and point to an intermediate inclination of ∼40° and a PA^kin ∼ 140°, both consistent with the results presented in Tan et al. (2019) and derived from ALMA observations of the CO(1-0) molecular gas emission. 3D-Barolo fit results are reported in Fig. C.1. Significant residuals are observed in the innermost nuclear regions, due to the AGN driven outflow, and close to the spiral arms, probably due to the fact that 3D-Barolo uses a concentric rings structure instead of spiral models to reproduce the observed kinematics.

I07251 This target shows two interacting nuclei at a projected distance of ∼1.8 kpc, and very extended tails and streaming gas with velocities from ∼ − 150 km s⁻¹ (north-east) to ∼ + 150 km s⁻¹ (north-west). In the innermost nuclear regions, it presents disk-like motions, but with a kinematic centre not coincident with none of its two nuclei (Fig. A.5). I07251 flux distribution is very irregular, and cannot be modelled with elliptical isophotes.

We fitted the narrow Hα data cube with 3D-Barolo, following the I method strategy, fitting v_rot, σ, ϕ, and the kinematic centre position. Unfortunately, the data quality does not allow us to constrain the disk inclination in this target; we therefore performed the second step of the 3D-Barolo fit assuming a mean inclination of 52° (Bellocchi et al. 2013) as initial guess. The 3D-Barolo best-fit results are shown in Fig. C.2. Also for this target, we observe significant residuals in the velocity and velocity dispersion maps, due to the complex nature of this interacting system.

Fig. C.2. I07251 velocity maps and the 3D-Barolo disk kinematic best fit. In this target, the kinematic centre is shown with a black cross and the two ULIRG nuclei with red crosses. See Fig. C.1 for further details.

I10190 W I10190 shows two interacting nuclei at a projected distance of 7.2 kpc, and very extended tails and streaming gas with velocities from ∼ − 200 km s⁻¹ (south-east) to ∼ + 300 km s⁻¹ (north-west). Each nucleus shows disk-like kinematics on kiloparsec scales (Fig. A.7); however, both the stellar and gas kinematics in the vicinity of the E nucleus are strongly affected by the presence of the W system motions (see e.g. Fig. B). This limits the possibility to study the E nucleus kinematics.

We fitted the I10190 W gas kinematics with both I and II methods, obtaining totally consistent results (Fig. C.3). Also for this target, we found a slightly offset between the near-IR nucleus and the fitted kinematic centre (∼0.3″); significant residuals are also observed in the 3D-Barolo maps.

Fig. C.3. I10190 W velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for further details.

I12072 N This source presents two nuclei, at a projected distance of 2.3 kpc, and a very extended plume in the north-east (Fig. A.9). Disk-like motions are presumably associated with the N nucleus. The flux distribution is complex and does not allow us to perform a robust Isophote modelling. We therefore fitted the narrow Hα data cube with the 3D-Barolo, applying the I method. The fit results are reported in Fig. C.4. Also in this case, velocity and velocity dispersion maps present significant residuals, due to the complex nature of this ULIRG.

Fig. C.4. I12072 velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for further details.

I13120 The 3D-Barolo analysis is extensively presented in Sect. 3.5.

I14348 NE I14348 is an interacting ULIRG with two nuclei at a projected distance of 5.3 kpc. The NE nucleus presents disk-like motions, while the SW kinematics are dominated by a strong outflow pointing to the south-west; streaming motions along the north-east south-west direction are also present on scales of 10s kpc (Fig. A.11). We fitted with 3D-Barolo the kinematics in the NE nuclear region, applying the I method to infer, as a first step, the inclination of the disk. In fact, the complex nature of this interacting system prevents robust Isophote analysis, and no morphological information is available. We decided to limit the 3D-Barolo fit analysis to the innermost nuclear regions (∼3″ × 3″), in order to exclude the contribution from the large-scale streaming motions. As a consequence, the disk inclination we determined with the I method, 52° ±3°, has to be taken with caution. The fit results are reported in Fig. C.5. Significant residuals are observed in the velocity and velocity dispersion maps. The kinematic centre has been fixed at the position of the nucleus during the 3D-Barolo analysis, as no significant variations in fit results are observed considering the kinematic centre position as free parameter.

Fig. C.5. I14348 NE velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for further details.

I17208 This source presents disk-like kinematics and extended tidal tails with velocities from ∼ − 150 km s⁻¹ (north) to ∼ + 150 km s⁻¹ (south-west; see Fig. A.14). For this target, we decided to exclude from the 3D-Barolo analysis the more external regions, associated with streaming motions; as a result, I method does not provide robust constraints for the disk inclination. We therefore applied the II method, assuming a disk inclination of 40° (from Isophote analysis) as initial guess. The 3D-Barolo fit results are reported in Fig. C.6. The kinematic centre has been fixed at the position of the nucleus during the 3D-Barolo analysis; no significant variations in fit results are obtained adding the kinematic centre position as free parameter.

Fig. C.6. I17208 velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for details.

I20100 SE The I20100 ULIRG is an interacting system with two nuclei at a projected distance of 6.5 kpc. The SE nucleus presents disk-like motions, while the NW kinematics are more complex and probably affected by the streaming motions (Fig. A.17). We therefore decided to analyse with 3D-Barolo the kinematics in the vicinity of the SE nucleus only.

The results obtained with methods I and II are broadly consistent (see Table 2), and point to an inclination of ∼58° and a PA^kin ∼ 287°. The significant difference between gas and stellar PA^kin might suggest more complex kinematics in the innermost nuclear regions of this source; indeed, we observe significant residuals in the velocity and velocity dispersion maps (Fig. C.7) along the gas kinematic minor axis.

Fig. C.7. I20100 SE velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for details.

Appendix D: Extended sample of SB disk galaxies

In this section we introduce the additional individual ionised gas measurements of SB disk galaxies reported in Fig. 7. In particular, we selected the following measurements:

i) the Hα (narrow component) velocity dispersion of 33 LIRGs and one ULIRG with disk kinematics from Bellocchi et al. (2013, B13 in the figure), observed with VLT/VIMOS, and the Brγ velocity dispersion of seven LIRGs from Crespo et al. (2021, C21 in the figure), observed with VLT/SINFONI. All these systems are characterised by mean σ of the order of several 10s km s⁻¹; ii) the Hα velocity dispersion of 22 rotationally supported dusty galaxies of the cluster Cl0024+17 at z ∼ 0.4 (Johnson et al. 2016, J16 in the figure), observed with VLT/FLAMES. These systems, with their median v_rot/σ = 5 ± 2 and their SFRs likely enhanced by the effects of ram pressure, also tend towards higher values compared to MS galaxies at the same redshift; iii) the Hα velocity dispersion of eight (U)LIRGs at z ∼ 0.2 − 0.4 from Pereira-Santaella et al. (2019, PS19 in the figure), observed with the optical integral field spectrograph SWIFT. Most of them are interacting systems (six out of eight) and have relatively small v_rot/σ ratios (from 0.4 to 3.2); iv) the Paα velocity dispersion of three SB disk galaxies at z ∼ 0.15 presented in Molina et al. (2020, M20 in the figure), and observed with SINFONI; v) the Hα velocity dispersion of seven disk galaxies at z ∼ 2 from the SINS/zC-SINF AO Survey, presented in Förster Schreiber et al. (2018, FS18 in the figure): ZC400528, ZC406690, ZC407302, ZC410123, ZC411737, ZC413507, and ZC415876. These galaxies display an offsets by a factor of ≳4 in SFR from the MS in their Fig. 6, and σ₀ ∼ 30 − 60 km s⁻¹; vi) the Hα velocity dispersion of a SB disk galaxy at z = 2.028, SMM J0217-0503b, merging with an AGN source (at a projected distance of ∼15 kpc), from Alaghband-Zadeh et al. (2012, A12 in the figure), and observed with VLT/SINFONI; vii) the Hα velocity dispersion of a SB disk galaxy at z = 2.24, SHiZELS-14, from Cochrane et al. (2021, Co21 in the figure). This source has been observed with the adaptive optics assisted VLT/SINFONI spectrograph; viii) the [OIII] (narrow component) velocity dispersion of a sub-millimetre bright ‘disky’ galaxy hosting a broad line quasar, SMM J1237+6203 at z = 2.075, from Harrison et al. (2012), and observed with Gemini-North NIFS.

All these measurements have been obtained from integral field spectrograph data. For the B13, C21, J16, A12, and H12 sources, σ₀ was computed as a mean velocity dispersion across the galaxy extension in the FoV; PS19 σ₀ were inferred from GalPak^3D modelisation (correcting for the beam-smearing and line-spread-function broadening), which assumes a spatially constant velocity dispersion in the disk; FS18 σ₀ were measured along the kinematic major axis at the largest radii possible, away from the central peak caused by the steep inner disk velocity gradient; M20 and Co21 σ₀ were instead derived as mean velocity dispersion in the external regions of the sources, excluding the innermost regions affected by beam smearing.

All Tables

All Figures

Fig. 1. Examples of our multi-component Gaussian fit decomposition for two continuum-subtracted spectra of I13120, extracted from single spaxels at ∼1″ (top panels) and 2″ (bottom) south-east of the nucleus. The red curves show the best-fit solutions in the Hβ–[OIII] (left) and Hα–[NII] (right) regions; all emission lines are fitted simultaneously. The vertical solid lines mark the rest-frame wavelength of the emission lines, derived from the stellar velocities of the nuclear spectrum; the vertical dot-dashed blue lines mark instead the local stellar systemic (i.e. at the position of the spaxel from which the spectrum is extracted). These examples demonstrate the diversity of emission-line profiles observed in the FoV of a single target: the spectrum in panels a is dominated by broad and blueshifted Gaussian components, while the one in panels b is dominated by bright narrow Gaussians close to the systemic velocity (for all but the [OIII] lines).

In the text

Fig. 2. I13120 maps. Top: total Hα integrated flux (left), Hα centroid (v50, centre), and line width (W80, right) obtained from the multi-component Gaussian fit. Bottom: similar panels for [NII]; in the W80 map, the red lines indicate the wide biconical outflow along the north-south direction. Masked regions mark the spaxels contaminated by the presence of background and foreground sources and excluded as disturbing the data analysis. The first solid contour is 3σ, and the jump is 0.5 dex. The cross marks the nucleus. North is up, and west is to the right.

In the text

Fig. 3. I13120 velocity diagram and maps. Left panel: I13120 velocity shift ΔVj − FWHMj diagram (coloured by density in log-scale) for the individual Gaussian components used to model the emission line profiles in the data cube. The dashed red lines isolate the Gaussian components used to reconstruct the narrow Hα data cube, with |ΔVj|< 250 km s−1 and FWHMj < 400 km s−1. Right panels: narrow Hα flux distribution, velocity (v50), and line width (W80) maps, reported in the second, third, and fourth panels, respectively (see Fig. 2 for details).

In the text

	Fig. 4. I13120 PV diagrams along the galaxy major axis. Left two panels: PV maps of the total Hα and the narrow Hα emission. Right panels: extracted line velocity centroids and line width of the narrow Hα, as well as the stellar velocities along the same axis. 1″ corresponds to ∼0.6 kpc at the distance of I13120, as labelled in the third panel.
In the text

Fig. 5. PA diagrams. Left: comparison between morphological and kinematic major axis PAs for the stellar (yellow) and gas (red) components. The dashed blue line indicates the 1:1 relation; the blue-shaded region (±22°) includes 90% of the morpho-kinematic PA misalignments of a sample of 80 non-interacting galaxies from the CALIFA survey (Barrera-Ballesteros et al. 2015). Right: comparison between stellar and gas kinematic major axis PAs. The dashed line indicates the 1:1 relation; the shaded region (±15°) includes 90% of the kinematic PA misalignments of non-interacting CALIFA galaxies (Barrera-Ballesteros et al. 2015).

In the text

Fig. 6. I13120 narrow Hα 3D-Barolo disk kinematic best fit of the moment 0, 1, and 2 (first to third rows) and PV diagrams along the minor and major disk kinematic axes (bottom). The black curves in the velocity dispersion map obtained from the data (third row, left panel) identify the region from which σ0 is extracted. In the PV diagram along the major axis (bottom right) and minor axis (left), data are indicated with a grey-scale map and blue contours, while best-fit models are shown with red contours.

In the text

Fig. 7. Gas velocity dispersion σ0 evolutionary trend and correlations. Top: velocity dispersion, σ0, as a function of z for the PUMA sub-sample (red dots) and other individual ionised gas measurements from the literature, distinguishing between the different samples presented in Appendix D, as labelled. The Übler et al. (2019) evolutionary trend of MS galaxies is shown with a solid curve (shaded area: 1σ scatter around the average trend). Middle: σ0 as a function of the SFR for all targets already reported in the first panel and the KMOS3D galaxies (grey points; Übler et al. 2019). K18 models and a linear fit (with scatter) are also reported, as labelled. Bottom: σ0 normalised to the evolutionary trend of MS galaxies as a function of δMS for all targets already reported above. A linear fit is also reported.

In the text

	Fig. A.1. I00188 maps: Hα integrated flux (left), Hα centroid (v50, centre), and line width (W80, right) obtained from the multi-component Gaussian fit. The first solid contour is 3σ, and the jump is 0.5 dex. The cross marks the nucleus. North is up, and west is to the right.
In the text

	Fig. A.2. IZw1 maps. See Fig. A.1 for details.
In the text

	Fig. A.3. I01572 maps. See Fig. A.1 for details.
In the text

	Fig. A.4. I0589 maps. See Fig. A.1 for details.
In the text

	Fig. A.5. I07251 maps. See Fig. A.1 for details.
In the text

	Fig. A.6. I09022 maps. See Fig. A.1 for details.
In the text

	Fig. A.7. I10190 maps. See Fig. A.1 for details.
In the text

	Fig. A.8. I11095 maps. See Fig. A.1 for details.
In the text

	Fig. A.9. I12072 maps. See Fig. A.1 for details.
In the text

	Fig. A.10. I13451 maps. See Fig. A.1 for details.
In the text

	Fig. A.11. I14348 maps. See Fig. A.1 for details.
In the text

	Fig. A.12. I14378 maps. See Fig. A.1 for details.
In the text

	Fig. A.13. I16090 maps. See Fig. A.1 for details.
In the text

	Fig. A.14. I17208 maps. See Fig. A.1 for details.
In the text

	Fig. A.15. I19297 maps. See Fig. A.1 for details.
In the text

	Fig. A.16. I19542 maps. See Fig. A.1 for details.
In the text

	Fig. A.17. I20087 maps. See Fig. A.1 for details.
In the text

	Fig. A.18. I20100 maps. See Fig. A.1 for details.
In the text

	Fig. A.19. I22491 maps. See Fig. A.1 for details.
In the text

Fig. B.1. Position-velocity diagrams along the galaxy major axis for all PUMA systems that show i) a well-defined velocity gradient along the major axis and ii) a peak in the velocity dispersion diagram close to the position of the nucleus. These two conditions provide initial evidence for rotation-dominated kinematics. All systems whose gas (stellar) kinematics satisfy the two conditions are marked with a blue (red) check mark in the top-right corner of the velocity dispersion panel; on the contrary, the systems that do not satisfy at least one of the two conditions are marked with a cross symbol; more uncertain kinematics are marked with a question mark. The systems for which neither the gas nor the stellar components satisfy the two conditions are not reported in this figure. Line velocity centroids and line widths of narrow Hα, total Hα, and stars are reported for each target, as labelled. See Fig. 4 for further details about the individual panels.

In the text

	Fig. B.1. continued.
In the text

Fig. C.1. IZw1 velocity maps and the 3D-Barolo disk kinematic best fit. Top panels: Narrow Hα (left) and pPXF stellar (right) velocity maps. The magenta line identifies the major axis PA measurement, computed within the black box region; the black cross identifies the nucleus. Second to fourth row panels: Comparison of 3D-Barolo data and model moment maps, as labelled in the figure. In the intensity and velocity maps, we report the major axis PA (dashed line) and the position of the kinematic centre (black cross); the green curves in the velocity maps represent the zero-velocity axis. In the dispersion map on the left, the black curves show the region from which the median σ0 value has been derived. Bottom panels: PV diagram along the major axis (right) and minor axis (left) for both the data (grey map and blue contours) and the best-fit model (red contours; the yellow dots associated with individual concentric rings are used to model the data).

In the text

	Fig. C.2. I07251 velocity maps and the 3D-Barolo disk kinematic best fit. In this target, the kinematic centre is shown with a black cross and the two ULIRG nuclei with red crosses. See Fig. C.1 for further details.
In the text

	Fig. C.3. I10190 W velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for further details.
In the text

	Fig. C.4. I12072 velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for further details.
In the text

	Fig. C.5. I14348 NE velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for further details.
In the text

	Fig. C.6. I17208 velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for details.
In the text

	Fig. C.7. I20100 SE velocity maps and the 3D-Barolo disk kinematic best fit. See Fig. C.1 for details.
In the text