© A. Calabrò et al. 2019

1. Introduction

Starburst galaxies (SBs) host the most powerful star-formation events in the Universe, with star-formation rates (SFR) that, following to the definition of Rodighiero et al. (2011), are at least four times higher compared to the average galaxy population at a given redshift, identified by the star-forming main sequence (MS; Brinchmann et al. 2004; Noeske et al. 2007; Elbaz et al. 2007; Daddi et al. 2007a). In the local Universe these systems are called ultra-luminous infrared galaxies (ULIRGs) and have a global infrared luminosity higher than 10¹² L_⊙, due to their high stellar production. ULIRGs are characterized by peculiar and rather homogeneous properties: they show very compact starbursting cores of dense gas, with typical diameter sizes of 0.5−2 kpc, and may contain double nuclei (Genzel et al. 1998; Díaz-Santos et al. 2010). As a consequence of the enhanced SFR density, these cores are also highly dust-rich and completely obscured in the optical/UV, with A_V ranging 5−50 mag in a foreground dust screen model (Genzel et al. 1998; Goldader et al. 2002).

SBs represent a small fraction of the whole star-forming (SF) galaxy population (2−3%), almost constant with redshift, according to Sargent et al. (2012, 2014) and Schreiber et al. (2015), but they constitute a key event in the life of a galaxy. Indeed, at z ∼ 0 the only viable option to explain these concentrations of gas, dust and SFR is through highly dissipative major merger events (Sanders & Mirabel 1996), in which the gas of colliding galaxies loses angular momentum and energy during the interaction, falling into the coalescing center of the system. Here it serves as fuel for the starburst and for the growth of a supermassive black hole in the center. This scenario is sustained by studies of both local ULIRGs and higher redshift merging and starburst galaxies, showing the presence of an AGN with optical (Ellison et al. 2011), mid-IR (Daddi et al. 2007b; Weston et al. 2017; Satyapal et al. 2017; Goulding et al. 2018), X-ray (Brusa et al. 2010; Aird et al. 2019) and radio (Best & Heckman 2012; Chiaberge et al. 2015) diagnostics. The starburst activity and subsequent AGN feedback can cause gas consumption and removal through powerful winds (Sanders et al. 1988; Silk & Rees 1998), leading to a passively evolving elliptical galaxy (Kormendy & Sanders 1992; Springel et al. 2005a; Hopkins et al. 2008a).

Hydrodynamical simulations indicate that the morphology of merger remnants is consistent with early type galaxies (ETGs), suggesting the scenario that gas-rich major mergers are primarily responsible for the formation of the massive ellipticals (Barnes & Hernquist 1991; Mihos & Hernquist 1994, 1996). This evolutionary sequence is strengthened by observations of low surface brightness tidal tales and residual interacting features in local ETGs (Duc et al. 2013), which may come from major mergers. Other studies have shown that ULIRGs lie on or close to the fundamental plane (effective radius-velocity dispersion-surface brightness) of intermediate-mass ellipticals, S0 galaxies and bulges of local spirals, suggesting that these systems are the final step of gas-rich merger episodes (Genzel et al. 2001).

The morphological transformation during the merger happens in different phases that are well represented by the so-called Toomre sequence (Toomre & Toomre 1972). Following the result of their simulations, several studies have tried to characterize this sequence by looking at different physical properties of the nuclear regions of merger galaxies, but found no correlations with evolutionary stages (e.g., Laine et al. 2003), apart from noticing that latest interaction times, during or after the coalescence, have among the highest infrared luminosities (L_IR) (Laine et al. 2003; Haan et al. 2013). The absence of clear correlations in these works can be interpreted either with the difficulty of identifying correctly the merger phase from the optical morphology (which becomes even more problematic beyond the local Universe) or that other parameters are driving and tracing this transformation. For example, Gao & Solomon (1999) and Leech et al. (2010) studied a sample of local (U)LIRGs with double merging nuclei. They found lower molecular gas masses and higher star-formation efficiency and gas excitation (probed by the CO(3–2)/CO(1–0) line ratio) with decreasing separation between merger components, that is, toward more advanced merger stages. However, since their sample only includes interacting pairs more distant than ∼1 kpc, this result is limited to relatively early stage mergers and does not necessarily apply also to coalesced systems. Furthermore, some spatially resolved studies on local (U)LIRGs have found that there is a higher contribution of shocks accompanied by an increased velocity dispersion of the gas toward later merger stages (Monreal-Ibero et al. 2006, 2010; Rich et al. 2014, 2015).

While in the local Universe the relation between starbursts and mergers is well settled observationally (e.g., Murphy et al. 1996; Luo et al. 2014), at higher redshift the interpretation becomes less clear. Due to the enhanced gas fractions and disk instabilities of high redshift star-forming galaxies, mergers as a result might increase SFRs less dramatically (Fensch et al. 2017), and starburst galaxies may be triggered by anomalous gas accretion events (Scoville et al. 2016). However, other studies have shown that the most extreme starbursts are still merger-driven, displaying disturbed morphologies (Elbaz & Cesarsky 2003) and increased star-formation efficiencies compared to MS galaxies (Sargent et al. 2014; Silverman et al. 2018a,b, 2015). Similarly, the connection between mergers and AGNs is still debated. Even though it may still hold for the most luminous cases (Combes 2003), some studies (focused especially on optical wavelengths) do not find systematic differences in merger fraction and galaxy distortions between active and non-active systems (Cisternas et al. 2011; Kocevski et al. 2012), and not all AGNs are triggered by mergers, according to, for example, Draper & Ballantyne (2012) and Ellison et al. (2015).

The solution to these long-standing problems is further complicated due to the intrinsic faintness of interacting features (e.g., tidal tales, bridges) and their elevated obscurations, which hamper their identification and physical characterization. Puglisi et al. (2017) have selected a sample of starbursts (SFR  ×  4 above the MS) at z ∼ 1.6, showing that optical lines, including Hα and Hβ, only probe a small nearly-unattenuated component of the galaxies, accounting for ∼10% of the total SFR_IR. Therefore, a possible solution to study the properties of the dusty starburst population at high redshift is to observe at longer wavelengths, targeting near-infrared rest-frame lines, such as the Paschen line series of hydrogen. Current ground-based spectrographs can observe Paβ (the second brightest Paschen recombination line) in K band up to a redshift of ∼0.9. Motivated by this idea, we followed-up of a sample of 25 SBs at 0.5 <  z <  0.9 with FIRE (folded port infrared echellette), a near-IR spectrograph mounted at the Magellan telescope, in order to detect optical and near-IR lines ranging 0.6 <  λ <  1.3 Å rest-frame. This sample was presented in Calabrò et al. (2018), which hereafter is referred as Paper I.

By comparing Hα, Paβ fluxes and bolometric IR luminosities, we found that the attenuation properties of intermediate redshift SBs are not consistent with the predictions of dust-foreground extinction models, but rather with those of a geometrical model in which dust and stars are homogeneously mixed (Paper I). We also found that they are highly obscured on average, with median A_V, tot = 9 mag in a mixed model, while independently derived dust-column densities suggest A_V, tot even higher for a large fraction of them, up to 75 mag. This means that they have extremely obscured cores and that optical and near-IR lines only probe an external skin containing a fraction (∼10−30%) of the total SFR. Even more, we argued that the presence of optically thick cores are themselves striking evidence of the merger origin of intermediate-z SBs (like in local ULIRGs), as no other mechanisms are known to produce such large obscuring column densities of gas and dust. This can be used to identify other mergers in the high-redshift Universe (Calabrò et al. 2018), where other methods based on morphology become unfeasible.

Despite their higher average obscurations compared to MS galaxies, we found that our sample, while showing a nearly constant Paβ/Hα ratios, spans a large range of A_V, tot between 2 and 30 (see Fig. 4 in Paper I), thus forming a sequence of increasing attenuations. This indicates that a substantial intrinsic diversity should exist among them, possibly related to different phases of the merger, to the gas properties (i.e., total gas mass, gas fraction) or to the morphology of preexisting colliding galaxies. In alternative, the sequence may be driven also by the orbital geometry of the merger. For example, Di Matteo et al. (2008) have shown that the inclination of the two encounters affects both the duration and intensity of the star formation, regulating the amount of gas that is funneled toward the coalescing center.

In this paper we investigate the origin of the obscuration sequence that was found in Paper I, analyzing the attenuations A_V, tot of our starbursts in relation with other physical properties derived from our Magellan spectra and from available multiwavelength images. In Sect. 2 we describe how the starbursts are selected from the parent photometric catalog along with our Magellan-FIRE observations, spectra reduction (which includes also the analysis of public optical spectra available for a fraction galaxies) and subsequent emission line fluxes (and equivalent width) measurements. This will be followed by a description of the morphological properties, radio size measurements and AGN identification within our sample. In Sect. 3 we will show the main results of this work, which are discussed extensively in Sect. 4. A summary with conclusions will be presented in the last Sect. 5. We adopt the Chabrier (2003) initial mass function, AB magnitudes, and standard cosmology (H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, Ω_Λ = 0.7). We also assume by convention a positive equivalent width (EW) for emission lines and a negative EW for lines in absorption.

2. Methodology

In this section we outline the general starburst selection procedure, and then describe the observations, the spectral reduction and the emission line measurements that were performed for a representative subset of 25 of them. We will then present their morphological classification, radio size measurements and AGN identification procedure.

As described in Paper I, we already calculated for this subset the absolute dust attenuation in V band toward the center (=A_V, tot) in two steps. First, we considered the ratio (SFR_IR/SFR_Paβ, obs), where SFR_IR is the intrinsic SFR from the infrared and SFR_Paβ, obs is derived from the observed Paβ luminosity, adopting an intrinsic ratio Paβ/Hα = 0.057 and a standard Kennicutt et al. (1994) calibration, valid for case B recombination and T_e = 10⁴ K. A_Paβ, IRX can be eventually inferred from that quantity as =2.5 × log₁₀ (1 + SFR_IR/SFR_Paβ, obs). In the final step, the attenuation in V band is obtained by solving for A_V, tot the following equation valid for a mixed geometry model of dust and ionizing stars:

$$\begin{array}{c}\hfill \frac{{\mathrm{SFR}}_{\mathrm{Pa}\beta ,\mathrm{obs}}}{{\mathrm{SFR}}_{\mathrm{IR}}}=\frac{{log}_{10}(e)}{0.4}\times \left(\frac{1-{10}^{-0.8A_{\mathrm{tot}}(\mathrm{Pa}\beta )}}{2A_{\mathrm{tot}}(\mathrm{Pa}\beta )}\right)\end{array}$$(1)

where A_tot(Paβ) is the total absolute attenuation at Paβ toward the center defined as k(Paβ) × A_V, tot/R_V. In the last expression, k(Paβ) corresponds to the local extinction at Paβ, assuming a Cardelli et al. (1989) law (R_V = 3.1). More details on the derivation and the implications of Eq. (1) can be found in Paper I.

2.1. Sample selection and Magellan FIRE observations

The starburst candidates were selected from the IR+(sub)mm catalog of Jin et al. (2018) by first requiring the spectroscopic redshift (coming from optical surveys, Salvato et al., in prep.) to be in the range between 0.5 and 0.9. This guarantees that Paβ and Hα fall, respectively, in the K and YJ bands, thus within the wavelength coverage of FIRE (0.82 μm–2.4 μm). For two galaxies, we found that their previous spectroscopic redshift measurements were not correct. Even though they are outside of our selection range, we did not remove them since they satisfy the other selection criteria.

Secondly, we required SFR_norm to be more than a factor of four higher than SFR_MS, z = 0.73, as in Rodighiero et al. (2011) ¹, where SFR_norm is the total SFR (SFR_tot = SFR_IR + SFR_{UV, unobscured}), normalized to the median redshift of the sample (0.73) using the evolution of Sargent et al. (2014):

$$\begin{array}{c}\hfill {\mathrm{SFR}}_{\mathrm{norm}}=({\mathrm{SFR}}_{\mathrm{IR}}+{\mathrm{SFR}}_{\mathrm{UV},\mathrm{unobs}})\times {\left(\frac{1+0.73}{1+z}\right)}^{2.8},\end{array}$$(2)

while SFR_{MS, z = 0.73} is the SFR of the main sequence at redshift 0.73 from Schreiber et al. (2015). This latter was found in Paper I to agree well with the MS derived independently for our parent sample through a running median on ten stellar mass bins from 10¹⁰ to 10¹² M_⊙. We chose this selection procedure in order to compare all the galaxies with a single MS (cf. Fig. 1 of Calabrò et al. 2018), instead of considering different relations for each galaxy redshift. However, the two procedures yielded the same number of starbursts among the whole sample. For the targets studied in this paper, the mean distance from the MS (dist_MS = log₁₀(SFR_tot/SFR_MS)) would change by 1.8%, thus it would not affect in any case the results.

The SFR_IR was derived following the methodology of Jin et al. (2018) and Liu et al. (2018), which is based on fitting the photometry from IRAC to radio 1.4 GHz with four components: a stellar Bruzual & Charlot (2003) SED (with age 200 Myr, constant star-formation history, solar metallicity Z_⊙, Chabrier IMF and Calzetti attenuation law), a warm+cold dust emission template from Draine & Li (2007) and Magdis et al. (2012), and a mid-infrared AGN SED from Mullaney et al. (2011) to separate the IR contributions of the AGN and SF. The unobscured UV-based SFRs (SFR_{UV, unobscured}) were calculated instead from the u-band magnitudes (Laigle et al. 2016), which probes the UV rest-frame at our redshift, following Heinis et al. (2014). Overall, the contributions of the AGN and of the UV-unobscured SFR to the total SFR are on average small (3% and 1%, respectively).

Finally, as a last requirement, we asked for M_* to be greater than 10¹⁰ M_⊙, where M_* comes from Laigle et al. (2016) and is computed at the photometric redshift. This condition ensures that, within the mass and redshift constraints adopted here, all SBs would be IR-detected with a S/N_IR >  5 (cf. Fig. 13 in Jin et al. 2018), thus we have a SFR complete sample of SBs. Then, for the whole sample, we considered the spectroscopic redshifts (when available) instead of the photometric values, however this does not significantly affect the stellar masses: the two estimates for the parent sample are in agreement within a 1σ scatter of 0.11 dex, compatible with the uncertainties reported by Laigle et al. (2016). Only for one SB analyzed in this work (ID 685067, z_spec = 0.37), the new stellar mass was remarkably lower (−0.56 dex), due to the large difference with previous photometric redshift (z_phot = 0.71). Therefore, dist_MS was even higher, strengthening its starburst selection.

This criteria yielded 152 SBs, 25 of which were observed during four nights at the Magellan 6.5 m Baade telescope (17−18/03/2017 and 22−23/03/2018). The observed 25 galaxies were chosen from the preselected SB sample according to a priority list. We preferentially observed sources close to bright stars (J <  19−20 mag), so as to facilitate target acquisition, although we eventually avoided blind offsets, since our galaxies are already sufficiently bright (peak magnitudes <19 mag) to be detected in ∼20−60 s in the good seeing conditions of those observations. In addition, we targeted intrinsically brighter sources first, maximizing $\text{SFR}/D_{\text{L}}^2(z)$ ratio (D_L is the luminosity distance), and assuming no prior knowledge about the dust attenuation of the system, which was set to 0 in all cases. This introduces a small bias in our selection toward the more massive objects. However, our galaxies span the full range of stellar masses above 10¹⁰ M_⊙. We refer to Fig. 1 in Paper I, where we presented the redshift, the SFR and the stellar mass distribution of our observed starbursts and our parent galaxy sample.

Our targets were observed with the single-slit echelle spectrograph FIRE (Simcoe et al. 2013), which has a wavelength range of 0.82−2.4 μm. We refer to Simcoe et al. (2013) for the full technical description of the instrument. We chose a slit width of 1″, (yielding a spectral resolution of R = 3000) to minimize slit losses (the average intrinsic FWHM angular size in K_s-band for our sample is ∼0.6″) and reduce the impact of OH sky emission. In all the cases, the slits were oriented along the semi-major axis of each galaxy, as determined from HST i-band images. Additionally, we benefited from good seeing conditions over all the four nights, with an average of 0.7″ and a minimum of 0.45″. The majority of our starbursts were observed in AB sequence, with single exposure times of 15−20 min. We chose single frame integrations of 20 min during the first run and 15 min in the second run, which significantly reduces saturation of OH lines in K band, thus helping spectral reduction in that band. We decided to double the integration times (completing the ABBA sequence) for galaxies with a lower S/N of the Hα+[N II] complex (based on real-time reduction), to improve the detection of fainter lines. In practice, doing an AB sequence is irrelevant for the spectral reduction, as the pipeline reduces each frame separately (see Sect. 2.2), though it allowed us to easily derive 2D emission line maps with the standard IRAF tasks IMARITH and IMCOMBINE.

2.2. Spectroscopic reduction

The spectra were fully reduced using the publicly available IDL-based FIREHOSE pipeline (Gagné et al. 2015). For each exposure, we used internal quartz lamps (one for each observing session) to trace the 21 orders of the echelle spectra and to apply the flat field correction. Then, the wavelength calibration was performed by fitting a low (one to five) order polynomial (depending on the spectral order) to ThAr lines of lamp exposures (taken close in time to the corresponding science frames). We checked that the residuals of the lines to the best-fit wavelength solution are less than 1 pixel in all the cases, and are <0.1 pixel for the majority of the orders. This translates into an average wavelength accuracy of Δλ/λ ≃  5  ×  10⁻⁵, nearly constant across the entire spectral range. Finally, the sky subtraction was applied independently for each single frame following the method of Kelson (2003). In this step, the OH lines in the spectra were used to refine the wavelength calibration.

The 1D spectra were extracted from the 2D frames using an optimal extraction method (Horne 1986). However, this procedure cannot be applied when there is a rapidly varying spatial profile of the object flux (Horne 1986), as in the presence of spatially extended and tilted emission lines. We used in these cases a boxcar extraction procedure, with a sufficiently large extraction aperture (always >1.3″) in order to include all the line emission from the 2D exposure. We adopted the boxcar extraction for three galaxies in our sample, which are the IDs 245158, 493881 and 470239. Given the agreement within the uncertainties between the fluxes measured with the two approaches for the remaining galaxies, the use of the boxcar procedure did not appear to introduce a systematic flux bias.

After the spectral extraction, we applied the flux calibration to each 1D extracted spectrum, using telluric spectra derived from the observations of A0V stars. Before dividing the object and telluric spectra in the pipeline, we could interactively refine the wavelength matching between the two by minimizing the rms of the product. However, slightly different times or airmasses between science and standard star observations at infrared wavelengths can produce non-negligible telluric line residuals, affecting the subsequent analysis. We found that this problem was more relevant in K-band, where strong telluric features are present in the observed wavelength ranges 19 950−20 250 Å and 20 450−20 800 Å. The residuals in these regions could produce artificial variations of the real flux up to a factor of 2, while it was less significant at shorter wavelengths (Y to H). In order to remove these artifacts, we followed the procedure described in Mannucci et al. (2001): we first considered a standard star at an airmass of ∼1.5 and calibrated it with two different stars observed at significantly lower and higher airmasses (e.g., 1.2 and 1.9). Then the two obtained spectra were divided, yielding a global correction function (which is different from 1 only in the regions of strong telluric features defined above) that applies to all the single-exposure spectra, each of them with a different multiplicative factor until the telluric line residuals disappear We fit a linear function in nearby regions free of telluric regions and emission lines, and then determined the correction function through minimizing the rms of the difference between the corrected spectrum and the aforementioned continuum fit. Finally, we combined (with a weighted-average) all the 1D calibrated spectra of the same object.

The error on the flux density f_λ obtained from the FIRE pipeline was checked over all the spectral range, analyzing the continuum of each galaxy in spectral windows of 200 Å and steps of 100 Å, masking emission lines. In each window, we rescaled the rms noise so as to have the ${\chi }_{\text{reduced}}^2=1$ when fitting the continuum with a low-order (≲1) polynomial. A spline of order 1 spanning the whole wavelength range was used as a correction function. This criterion, equivalently, ensured that the noise level matched the 1σ dispersion of the object spectrum in each window. Typical corrections were within a factor of 2, variable across the spectral bands.

Due to slit losses, variable seeing conditions and the spatial extension of our objects, which are typically larger than the slit width (1″), part of the total flux of the galaxies is lost. In order to recover the total absolute flux, we matched the whole spectrum to the photometric SED. This was done by applying a 5σ clipping and error-weighted average to the Magellan spectrum inside z++, Y, J, H and K_s photometric bands, and comparing the obtained mean f_λ in each filter to the corresponding broad-band photometry (Laigle et al. 2016). Since the SED shapes derived from the spectra are generally in agreement with the photometric SED shapes, we rescaled our spectra with a constant factor, determined through a least squares minimization procedure. The aperture correction factors for our sample range between 0.8 and 3, with a median of 1.4. They are subject to multiple contributions, such as slit position with respect to the object, seeing conditions during the target and the standard stars observations. The few cases in which the aperture correction was lower than 1 could be due indeed to a much better seeing of the standard star compared to the target observation. We remind that this procedure assumes that lines and continuum are equally extended, which is clearly an approximation. Spatially resolved near-IR line maps (e.g., with SINFONI) would be required to test possible different gradients of the two emission components, and to derive better total flux corrections.

2.3. Complementary optical spectra

A subset of our Magellan sample also has publicly available optical spectra: ten starbursts in our sample have been observed with the VIMOS spectrograph (Le Fèvre et al. 2003) by the zCOSMOS survey (Lilly et al. 2007), and their optical spectra are publicly available through the LAM website (cesam.lam.fr/zCOSMOS). They span the range 5550 <  λ <  9650 Å, which includes [O II]λ3727 Å, Hγ, Hβ and [O III]λ5007 Å lines for our galaxies. The spectral resolution is on average R = 600, constant across the whole range, while the noise level increases from the blue to the red part of the spectrum, due to the increased OH sky emission at longer wavelengths. Due to the absence of the noise spectrum in the public zCOSMOS release, we used instead a sky spectrum at the same resolution of VIMOS, rescaled with a spline to match the flux standard deviation in spectral regions free of emission lines, similarly to what has been done before on the Magellan spectra. Even though it is a first approximation, this procedure allows to reproduce the increasing noise in correspondence of OH lines, where strong sky-subtraction residuals are expected, and take into account the higher average noise level of the red part of the spectrum. For one galaxy in our sample (ID 493881), we took its optical spectrum from SDSS-III DR9 (Ahn et al. 2012), which spans a wider wavelength range 3600 <  λ <  10400 Å at higher resolution (R ∼ 2000), thus allowing a better sky subtraction and a higher S/N. In both cases, the spectra were already fully reduced, wavelength and flux calibrated, as described in the respective papers. We apply only an aperture correction by matching the observed spectrum to the photometric SED (Laigle et al. 2016) with the same methodology adopted for the Magellan spectra. However, we warn that there could be some mis-matches compared to our Magellan observations in the slit centering and orientation, as also in the seeing conditions, thus the spectra may not be exactly representative of the same regions.

2.4. Line measurements

We measured emission line fluxes, line widths and uncertainties on fully calibrated and aperture corrected spectra using the python anaconda distribution (Mark Rivers, 2002²) of the IDL routine MPFIT (Markwardt 2009). Given the FWHM resolution of FIRE for 1″ slit width (=100 km s⁻¹), all our emission lines are resolved, owing to intrinsically higher velocity widths.

We fit the lines with a single Gaussian on top of an underlying order-1 polynomial continuum. In each case, we required a statistical significance of the fit of 95%, as determined from the residual χ². When a single Gaussian did not satisfy the above condition, we used a double Gaussian (with varying flux ratio and same FWHM, in km/s, for the two components), which instead provided a better fit, placing its χ² within the asked confidence level. The flux uncertainties were derived by MPFIT itself, and they were always well behaved, with best-fit ${\chi }_{\text{reduced}}^2\simeq 1$. For non detected lines (i.e., S/N <  2 in our case), we adopted a 2σ upper limit. We remark that we were guided by the wavelength position, line width and flux ratio (for double line components) of Hα, which was always detected at >4σ. The Gaussian amplitude remained thus the only variable to constrain the fit for the other lines. However, we highlight that our detected emission lines have always high S/N ratios: Hα, [N II]λ6584 Å and Paβ were identified on average at 9.3, 8.4 and 7.4σ, respectively (lowest S/N were 4, 5.3 and 3.3 for the same lines).

We fit a double Gaussian for 12 galaxies in our sample. As we will see later in Sect. 4.1 by combining the informations of their 1D and 2D spectra, in six of them we interpret the two Gaussians as coming from different merger components. For the remaining galaxies, in two cases the lines are consistent with global rotation, while for the last four we were not able to derive firm conclusions, even though we favor the contribution of multiple system parts to their emission. In Appendix C, we show the 1D emission line fits for all our 25 starbursts, and we discuss in more detail the origin of double Gaussians line components.

We applied the stellar absorption correction on Balmer and Paschen emission lines, rescaling upwards their fluxes. In order to determine the level of absorption for these lines, we adopted Bruzual & Charlot (2003) synthetic spectra with solar metallicity and constant star-formation history for 200−300 Myr, which are the typical merger-triggered starburst timescales (Di Matteo et al. 2008). The current starburst activity imposed by our selection suggests that the final coalescence of the major merger occurred relatively recently, certainly within the last 200 Myr. Numerical simulations of major mergers with different masses and dynamical times indicate indeed that star formation stops within 100–200 Myr after the coalescence, even without AGN quenching (Springel et al. 2005b; Bournaud et al. 2011; Powell et al. 2013; Moreno et al. 2015). Averaging the results over this interval, we applied an EW_abs of 5, 2.5, 2.5 and 2 Å for Hβ, Hα, Paγ and Paβ, respectively. The EW_abs of Hβ and Hα are consistent with those adopted in previous works (e.g., Valentino et al. 2015), while it was not possible to make comparisons for Paγ and Paβ. In the same order, we estimated for those lines an average absorption correction of 35%, 7%, 26%, 13% of the total flux. If we allow an uncertainty of ±1 Å on the EW_abs correction of either Paγ and Paβ, this will produce variations on the final fluxes that are 6% on median average, and thus it will not significantly affect our results.

All the lines in the Magellan spectra, either in emission or in absorption, were analyzed based on the following steps. Firstly, Hα and [N II]λλ6548, 6583, which are the lines with the highest S/N, were fit together assuming a common linear continuum and a fixed ratio of the [N II] doublet of 3.05 (Storey & Hummer 1988). From this fit we derived the redshift of the galaxy, the intrinsic FWHM of Hα (in terms of velocity), and the flux ratio of the two Hα components for double Gaussian fits. The instrinsic total line widths were obtained by subtracting in quadrature the FIRE resolution width (100 km s⁻¹) from the best-fit observed FWHM. For double Gaussians, the total FWHM was calculated adding the single FWHM and the separation between the two component peaks, as this quantity is more representative of the whole system.

Then, the three parameters defined above were fixed and used to fit all other emission lines, including those in the optical zCOSMOS and SDSS spectra, after rescaling the FWHM to account for the different spectral resolutions. For the galaxies with the highest S/N of Paβ (>8σ), we verified that even without fixing its FWHM a priori, the fit yields a line velocity width consistent within the errors with the value found from Hα, indicating that our assumption is generally valid. Given the rms wavelength calibration accuracy (see Sect. 2.2), we allowed the line central wavelength to vary in the fit by 3σ, corresponding to 1.5 Å at 10 000 Å, and 3 Å at 20 000 Å. For each measured flux, we also added in quadrature an error due to the uncertainty of the absolute flux calibration. This additional uncertainty ranges between 5% and 10%, and is determined as the maximum residual difference between the average fluxes estimated from the photometry and from the aperture corrected spectra, among all the bands ranging from z++ to K_s. Finally, the equivalent widths of the lines were derived following its definition (=∫(F(λ)−F_cont(λ))/F_cont(λ)), where F(λ) is the best-fit Gaussian flux distribution and F_cont is the fit underlying continuum. Since the fluxes of Hα and Paβ were presented in Paper I, here we show in Table A.1 the FWHM of the lines, the EWs of Hα, Hδ and Paβ, the fluxes of [O III]5007 and Hβ that have been used in the BPT diagram.

2.5. Ancillary data

Almost all of our starbursts (24) were observed by HST-ACS in the F814W filter (Koekemoer et al. 2007) at an angular resolution of 0.095″ (∼0.7 kpc at z = 0.7). The UltraVISTA survey (McCracken et al. 2012) observed our galaxies in YJHK bands at a spatial FWHM resolution of ∼0.75″, comparable to the average seeing during our Magellan observations. Two galaxies in the subset have higher resolution (0.19″) F160W HST images from the DASH program (Momcheva et al. 2016). Finally, all our galaxies are well detected in radio 3 GHz VLA images (Smolčić et al. 2017) (owing a similar a spatial resolution of ≃0.75″), with an average S/N of 18. In Calabrò et al. (2018) we showed the HST F814W, H-band UltraVISTA and radio 3 GHz VLA cutout images for only ten galaxies in our whole sample, for reason of space. Therefore, we include here in Fig. 1 the same types of images for the remaining sample of 14 starbursts, all of which have been observed during the second Magellan run (we remind that galaxy ID 578239 was not observed by HST, so we do not show it).

Fig. 1. For each of the galaxies observed in our second observing run at Magellan from left to right: HST-ACS F814W images (FWHMres = 0.095″), H-band UltraVISTA cutouts (same f.o.v. and FWHMres ∼ 0.75″) and 3 GHz radio images from VLA-COSMOS 3 GHz large project (Smolčić et al. 2017; FWHM ∼ 0.75″). Cutout images for the sample of the first observing run were presented in Paper I.

2.6. Morphological classification

Even though the light emission of dusty starbursts at optical and NIR wavelengths might be still severely affected by dust, the high resolutions offered by HST F814W images allows us to investigate the global structure of these systems. In Paper I we showed that the morphology of 18 galaxies has been already classified by (Kartaltepe et al. 2010; K10), revealing a merger origin for the majority of them. Adopting the same criteria of K10, we classified visually the remaining six galaxies (one has no HST coverage), but the results do not change significantly: 61% of our total sample are major mergers, as revealed by their highly disturbed morphology, tidal tales and bridges, 23% are classified as minor mergers from the presence of only slightly perturbed structures (e.g., warped disks, asymmetric spiral arms, etc.) without large companions, 11% are classified as spheroidal or S0 galaxies and the remaining 5% as spirals. The major merger subset is additionally divided in five smaller classes according to their merger state (I: first approach, II: first contact, III: pre-merger, IV: merger, V: old-merger or merger remnant), following K10.

However, we remind that the merger recognition and, even more, the merger stage classification from the optical morphology is very uncertain and more difficult at higher redshifts, due to lower resolution and to surface brightness dimming, which hampers the detection of faint tidal tales or interacting features, especially after the coalescence. The galaxy ID 245158 represents a show-case example of this uncertainty: it has been classified as spiral or minor merger from its global morphology, but it clearly shows a double nucleus in the central region, further confirmed by a double component Hα in the 2D and 1D spectrum, indicating rather an ongoing merger system.

2.7. Radio size measurements

The radio continuum emission has been used as a dust-free tracer of the SFR in galaxies in the absence of contamination from an AGN (e.g., Condon 1992). Since all our galaxies do not show either radio jets or radio flux excess compared to that expected from the SFR only (as we will show in Sect. 2.8), we used 3 GHz VLA images for measuring the FWHM size of their starburst cores, where the bulk of star-formation is taking place. For each SB, we used GALFIT (Peng et al. 2010) to fit a 2D function, convolved with the VLA synthesized beam, to their radio emission. We tested several 2D profiles, which include a Gaussian, a Sérsic function and the VLA beam itself, requiring a significance of the fit (probed by the χ²) of at least 95% confidence level, as for the emission line measurements. In addition to the χ² analysis, all the GALFIT residuals of the fit (original-model) were checked by eye inspection, and the excluded fittings have always worst residuals. The best-fit profiles obtained for our sample are summarized as follows.

Firstly, a single 2D Gaussian with varying FWHM, axis ratio and position angle provided the best fit for 13 starbursts. In one case (ID 470239), the required conditions were obtained only by fitting a single Sérsic profile with varying parameters, but its half light diameter (calculated as 2 × r_e, with r_e the effective radius) was only 3% different from the total FWHM of a single Gaussian fit, thus we assumed the latter as the final value. We also tried to fit a single and double Sérsic profile for all the other sources. However, given the larger number of parameters of this profile and the limited VLA resolution, we did not obtain convergence for the majority of them, or the resulting ${\chi }_{\text{reduced}}^2$ were too high.

A double 2D Gaussian was required by three galaxies (ID 245158, 412250 and 519651), allowing to resolve them and measure single components FWHM and their relative separation.

As a third option, fitting the VLA synthesized beam yielded the best solution for six galaxies, which are then unresolved with current resolution (0.75″).

Finally, a single 2D Gaussian with fixed axis ratio and position angle (to 1 and 0, respectively) were used for two sources (ID 578239 and 685067). We remark that, in case the 95% significance level of the fit was satisfied with either this or the previous approach (as in the case of some very compact sources), we adopted the Gaussian solution only if the associated χ² probability level was at least double compared to the fit with the VLA beam.

As shown later, for a few galaxies we measured angular sizes that are much smaller than the synthesized beam FWHM (∼0.75″), down to ∼0.2″ and to a physical scale of 1 kpc. To demonstrate that it is possible to reliably determine the sizes even for these extreme, compact sources, we show in Fig. 2a comparison between the GALFIT residuals obtained when fitting the image with a Gaussian (convolved to the VLA beam) (upper row) and with the radio synthesized beam itself (bottom row). It is evident that a Gaussian provides a better fit of the original source profile and a cleaner residual compared to the VLA beam alone.

Fig. 2. GALFIT fitting of the radio VLA image (3 GHz) of the galaxy ID 685067 (z = 0.37) with the VLA synthesized beam in the upper row and with a Gaussian profile (convolved with the beam) in the bottom row. In horizontal sequence are shown, from left to right, the original image, the fitted model and the residual (original-model). For this galaxy, we derive an angular size of 0.20 ± 0.04″ (the pixel scale is 0.2″ pixel−1), which corresponds to a physical size of 1.06 ± 0.19 kpc. We notice that here GALFIT converges when fixing the position angle (PA) and axis ratio (q) parameters, to 0 and 1, respectively (Table A.1). This example illustrates the possibility to reliably measure the radio sizes of our objects even when they are smaller than the FWHM resolution of VLA (0.75″). In this case, the difference between the two models is recognized by looking at the residual images (i.e., original-model).

The uncertainties on the sizes were recalculated for all the starbursts from their radio S/N, using the fact that better detected radio sources also have the smallest radio size uncertainties, as shown for example in Coogan et al. (2018). We used then the same formulation as:

$$\begin{array}{c}\hfill FWH{M}_{\mathrm{err}}\simeq 1\times \frac{FWH{M}_{\mathrm{beam}}}{S/N},\end{array}$$(3)

where FWHM_beam is the circularized FWHM of the VLA synthetized beam, and the multiplying coefficient was determined from simulations, following Coogan et al. (2018). All the size measurements with relative uncertainties and the method used for their determination are included in Table A.1.

Since our galaxies are well detected in radio band (average S/N of 18) and the VLA synthetized beam is well known, we always obtained a good fit for the resolved sources. Among them, we were able to fit a double Gaussian for three objects. In these cases, their total FWHM (adopted throughout the paper) were determined as the sum of the average single FWHM sizes and the separation between the two components. However, we will also consider the single sizes in some cases, such as in Fig. 10. This finding suggests that also some other galaxies may represent double nuclei that are blended in 0.75″ resolution VLA images. For the unresolved galaxies instead (i.e., those fitted with the VLA beam, as explained above), we adopted a 3σ upper limit on their FWHM. Within the most compact starbursts, some of them may be affected by pointlike emission from an AGN, which decreases artificially the observed size. However, we tend to discard this possibility since, as we will see later, none of our AGN candidates show a radio-excess compared to the radio emission due to their SFR.

2.8. AGN identification

We started searching for AGN components in the mid-IR. Through the multicomponent SED fitting of IR+(sub)mm photometry (described in Sect. 2.1) we detected at >3σ the dusty torus emission component for a subset of 12 SBs. The significance of the detection was derived from the ratio between the total best-fit dusty torus luminosity (=L_AGN, IR) and its 1σ uncertainty, inferred as the luminosity range (symmetrized) yielding a variation of the ${\chi }_{\text{red}}^2\le 1$ with respect to the minimum value of best-fit (Avni & Bahcall 1976). More details about the torus estimation method are described in Liu et al. (2018) and Jin et al. (2018). Among the 12 mid-IR AGNs, we detected the dusty torus emission at high significance level (>5σ) for six starbursts (ID 777034, 519651, 222723, 232171, 466112, 894779), while for the remaining objects we obtained a lower significance ranging 3σ <  L_AGN, IR <  5σ (see Table A.1). The SED fitting of all the galaxies can be found in the Appendix B.

Within the sample of IR-detected AGNs, six galaxies (ID 777034, 222723, 232171, 635862, 578239, 911723) were also detected in X-rays at more than 3σ by XMM-Newton, Chandra or NuStar (Cappelluti et al. 2009; Marchesi et al. 2016; Civano et al. 2015). Throughout the paper, we will consider the X-ray luminosities L_X measured by Lanzuisi et al. (2017), integrated over the energy range 2−10 keV. To estimate the contribution of star-formation to the total intrinsic L_X, we used the relation between SFR and L_X, SFR of Mineo et al. (2014), rescaled to a Chabrier IMF and applying a correction factor of 0.6761 to convert the X-ray luminosity from the 0.5−8 keV to the 2−10 keV band.

We remind that the column densities N_H inferred from their hardness ratios (Lanzuisi et al. 2017) are consistent with those derived from the dust attenuations (toward the centers) assuming a mixed model (Paper I), suggesting that also the X-ray emission is coming from the nucleus, where the AGN is expected to be located. Furthermore, all our starbursts do not show radio jets in VLA images, and do not have significant radio excess than expected from their SFR, assuming a typical IR-radio correlation with q_IR = 2.4, as in Ibar et al. (2008), Ivison et al. (2010) and Liu et al. (2018).

3. Results

The spectra that we obtained at the Magellan telescope, along with longer wavelength radio images, provide us key information to understand both the attenuation sequence and the variety of morphological classes of our starbursts. First of all, since dissipative mergers are able to funnel the gas from the large scales of Milky-way-like disks (∼10 kpc) to sizes that are more than one order of magnitude smaller (Di Matteo et al. 2005), it is useful to analyze the characteristic star-forming sizes of our starbursts. Besides this, from the galaxy integrated Magellan spectra we can study together the excitation and kinematic state of the gas, and the aging of the stellar population in the outer starburst cores, traced respectively by the [N II]/Hα ratios, the intrinsic (resolution corrected) line velocity widths of single 1D Gaussian components (=FWHM_line, which is a proxy for the velocity dispersion in the system) and the Balmer (or Paschen) line equivalent widths (=EW_Hα, Paβ).

In Fig. 3 we present the main result of this analysis, showing that the FWHM_radio, the N2 parameter, the FWHM_line and the EW_Paβ are all correlated to the total dust attenuation A_V, tot, which is used here as the reference quantity for comparison. This suggests that our starbursts can be described as, at first order, a one-parameter sequence: similar correlations at different significance levels are indeed found also when comparing on a single basis each pair of the above parameters.

Fig. 3. Left column: correlations of AV, tot with the radio size (top), and with the N2 index (bottom). We show with filled blue diamonds all the Magellan SBs that were used to derive the best-fit linear relation (blue continuous line) and the ±1σ dispersion (blue shaded area), while empty diamonds represent galaxies excluded from those calculations. The latter comprise the four outliers discussed in the text (ID 303305, 345018, 472775, 545185) and all the upper limits in the AV, tot vs. FWHMradio plot. In the corners, we show in black the equation of the linear fit (which includes 1σ error of the two best-fit parameters), and in gray the Spearman correlation coefficient (R) with the corresponding p-value (p), the reduced chi-square of the fit (${\chi }_{\text{red}}^2$), and the 1σ scatter of our SBs around the best-fit line, all of which do not take into account upper limits and the four outliers mentioned above. For comparison, the linear fit and 1σ dispersion including the four outlier galaxies are highlighted with a gray continuous dashed line and two dotted lines of the same color. Right column: correlations of AV, tot with the line velocity width (top), and with the EW(Paβ) (bottom). In the last diagram, four galaxies without EW(Paβ) measurements are not considered.

We tested these correlations using three different approaches with all the available data, excluding from the calculations only the upper limits and missing EW(Paβ) measurements. Firstly, we calculated the Spearman rank correlation coefficient R (the higher R, the stronger the correlation) and the corresponding p-value, which represents the probability of obtaining an equal (or stronger) R if no correlation exists. We defined a threshold probability of 5% to accept the correlation. Overall, we found that the p-values are nearly always lower than 0.05, meaning that the correlations are significant according to our criteria. In only one case (EW_Paβ vs. FWHM_radio) we determined a slightly higher p-value of 0.1 (thus a higher probability of no correlation), which could be partly affected by the lowest number of data (i.e., lowest statistics) available here compared to the other diagrams. However, the other methods indicate instead a stronger physical connection between the two quantities.

In the second approach, we fit the galaxies in each diagram with a linear relation (in log–log space, except for the last diagram where the y-axis is in linear scale), by using an orthogonal distance regression procedure (ODR), which allows to take into account measurement uncertainties in both axis (we discuss later possible outliers or different fitting functions). We determined the S/N of the angular coefficient (i.e., how much it differs from 0), finding significant correlations at more than 3σ in eight cases, while they are less strong (2 <  S/N <  3) for the remaining two diagrams. In the four correlations shown in Fig. 3, we obtained a significance of 5.8, 5, 4.3 and 3.65σ for A_V, tot vs. N2, FWHM_line and EW_Paβ, respectively. With this method, we also determined the 1σ dispersion of our data with respect to the best-fit linear relation.

Finally, we also performed Monte Carlo simulations: for each relation, we run 100 k simulations, removing each time at random 20% of the points, recalculating the significance of the correlation using our second approach. We then estimated the rate (∼probability) at which such correlations completely disappear with a significance falling below 2σ. This analysis allows to test the systematics and scatter of the correlations, ensuring they are robust and not driven by a few outliers. Overall, we found low probabilities (less than 5%) to obtain a less than 2σ significance when removing a random 20% of the galaxies, indicating that our correlations do not cancel out and are not found by chance. In the four diagrams of Fig. 3, we obtained probabilities of 0.028%, 0.001%, 4.7% and 0.7%, in the same order as above.

We remind that A_V, tot were determined from the Paβ observed fluxes and the bolometric L_IR, assuming a mixed model geometry. As explained in Paper I, for one galaxy in our sample where Paβ resides in an opaque atmospheric region (ID 245158), we estimated the attenuation through the Paγ line, adopting a flux ratio Paβ/Paγ = 2.2. This is the average expected observed ratio for all the attenuation values in our range, assuming either a mixed model or a foreground dust-screen geometry, and it was verified by nine starbursts with simultaneous detection of Paγ and Paβ.

However, for four galaxies in the sample (ID 303305, 500929, 893857 and 232171) we did not detect either Paβ or Paγ, thus in these cases we derived A_V, tot in a similar way from their Hα fluxes (so to avoid upper limits), adding a representative error of 0.1 dex determined from the remaining sample as the scatter of the correlation between Paβ and Hα based A_V, tot. We also verified that including the upper limits in the calculations does not significantly alter the fitted trends. Hereafter, we discuss in detail on a single basis the most important findings.

In the first (top-left) panel of Fig. 3, the FWHM radio sizes, while spanning a wide range from less than 600 pc to ∼12 kpc, are tightly anticorrelated to the dust obscuration level A_V, tot (R = −0.6, p-value = 0.007, and a scatter of 0.26 dex). Toward the smaller sizes and higher obscurations (A_V, tot >  20 mag), three galaxies are unresolved with VLA, thus they may be actually closer to the best-fit linear relation derived from the remaining sample. In this diagram, X-ray detected AGNs are found both at small and large radii, and have a similar distribution compared to the other galaxies, suggesting that radio size measurements and hence the result in Fig. 3 are not contaminated by AGNs.

In the last three panels of Fig. 3, the [N II]/Hα ratio, the line velocity width (FWHM_line) and the EW of Paβ are also correlated to the total attenuation at more than 3σ significance level (R coefficients and p-values are 0.51(0.009), 0.48(0.015) and −0.46(0.034), respectively).

We used the EW of Paβ for comparison since, being at longer wavelength, it is more representative of the whole system, allowing to probe a larger fraction of starburst cores if a mixed geometry holds. However, in the Appendix B we show that EW(Paβ) is tightly correlated to the EW of Hα, Hβ and Hδ, all of them being strongly sensitive to the age of the stellar population (at fixed SFH), thus similar results are obtained also if choosing a different line for the EW.

We also notice that four galaxies (ID = 303305, 472775, 345018 and 545185) are outside the 1σ dispersion of the best-fit relations in all the three diagrams (gray dashed and dotted thin lines). They show lower N2, FWHM_line, and higher EW than expected from their dust obscuration level. Alternatively, they have a larger A_V, tot for their N2, FWHM_line and EW values.

In order to understand the nature of these galaxies, we simulated 100 k different realizations of the last three diagrams of Fig. 3, with N2, FWHM_line, EW(Paβ) and A_V, tot of 25 galaxies distributed according to the best-fit relations and the corresponding 1σ dispersions. Then we computed the probability of having at least four galaxies (three for the last plot) with an orthogonal distance from the best-fit relation (gray dashed line) equal or greater than the four (or three) outliers described above. We found, in the same order presented above, a probability of 0.2%, 0.025% and 0.005%, indicating that those four galaxies are real outliers and cannot be simply explained by the 1σ scatter of the best-fit lines.

Given their deviant behavior, we excluded these outliers and derived again the best-fit relations, which are shown in Fig. 3 with a blue continuous line. We found on average a reduction of the 1σ dispersion (shown with a light blue shaded area) by ∼0.1 dex and a slight improvement of the correlation significance compared to the previous calculations. However, the best-fit linear equations are not significantly different, thus we give only the analytic expressions of this second fit where the outliers are not considered. The new results for the three diagrams, and all the diagnostics for the remaining seven correlations are presented in Table 1. We notice that the four divergent starbursts have an upper limit on their radio size, and are not outliers in other diagrams that do not involve A_V, tot, thus the latter are not affected by this analysis. A possible physical explanation of the diverging behavior of these four galaxies will be discussed in Sect. 4.4.

Finally, if we look at all the correlations in Table 1, we can notice the presence of a subset of quantities that correlate better than others. Apart from the previously discussed A_V, tot vs. FWHM_size, the line width, N2 and radio size are tightly and robustly correlated with each other. Indeed, from bootstrapping analysis, the probability that there is no correlation is less than 0.033%. As we will see in Sect. 4.2, this result hides a deeper physical link among them.

4. Discussion

The results presented in the previous section show that the wide range of attenuations measured in Paper I translate into a wide range of other physical properties, such as radio sizes, N2, velocity width, Balmer and Paschen EW, and even more, all these quantities appear to be connected to each other, defining a one-parameter sequence. In this section we propose a physical interpretation of this sequence, and show that the correlating properties considered before are consistent with being primarily reflecting different evolutionary merger stages. Then we discuss the role played by each parameter into this sequence.

4.1. Identification of early-phase, precoalescence mergers

A first guess for a physical understanding of what is guiding the large spread of properties comes from the morphology. Indeed we have already seen that our sample comprises mergers at different stages of evolution (MI to MV), though this classification is very uncertain and sometimes misleading, as shown in Sect. 2.6: the faintness of tidal tails and residual interacting features make systems at the coalescence difficult to recognize, while multiple optical components and double nuclei in HST images may just reflect the dust attenuation pattern rather than the true distribution of SFR and M_*. As shown in Sect. 2.7, a double Gaussian component fit on radio images allowed to resolve three sources, suggesting that they may be composed of two interacting nuclei. However, the limited resolution of VLA (0.75″ of FWHM) does not allow us to derive solid conclusions on the remaining sample, which might contain more close pairs. New maps and ALMA follow-ups would increase the resolution and hopefully resolve these blended precoalescence systems.

A complementary way to find close interacting pairs in early merger stages comes from the analysis of their 2D spectra. With that aim, we performed a crude sky-subtraction procedure: we subtracted the 2D spectra taken for the same object but at different positions along the slit (A and B, separated by 2.2″), in order to remove the sky lines and allow a visual inspection of the emission line profiles. For construction, the lines appear twice in each sky-subtracted frame and exactly with the same shape: one time with a positive flux (when the object is in position A), and the other with a negative flux (when the galaxy is in position B).

By looking at these line profiles, we identified interacting pairs by requiring one of the following conditions. Firstly, we looked for detached Hα line components along the spatial direction, coming from separated merger components located at different slit positions (e.g., ID 223715, 519651, 545185, 668738 in Fig. 4). Alternatively, we required tilted Hα line with two different inclination angles (based on visual inspection), indicating the presence of two emitting regions with independent kinematic properties, inconsistent with a single rotating disk (e.g., ID 245158, 493881, 470239 in Fig. 4). In our sample, we identified from the two above conditions seven close-pair precoalescence starbursts, which are shown in Fig. 4. For an additional source with a double radio emission component (ID 412250), one of the two nuclei was not falling inside the slit, thus it was not observable with FIRE. However, this SB should be considered a merging pair at the same level of the others.

Fig. 4. Close-up view of the Hα emission line profiles for the galaxies satisfying one of the two precoalescence criteria defined in Sect. 4.1. In each panel, as a consequence of the sky-subtraction procedure applied to the 2D spectra, the lines appear twice in different slit positions, in the first with a positive flux (in yellow) and in the second with a negative flux (in black). For each cutout, we show with two white arrows the slit position (s) and the dispersion direction (λ), which are slightly rotated due to the curved Magellan spectra. In order to clarify the classification criteria adopted in this work for finding precoalescence SBs, we highlight: (1) with continuous lines the different tilting angles of Hα line profiles (first three SBs); (2) with dotted ellipses the spatially separated Hα lines coming from different merger components (last four SBs).

As we can see in Fig. 3, the selected precoalescence starbursts are preferentially found at larger half-light radii, and all the systems with FWHM_radio >  6 kpc belong to this category. This result has two main implications. Firstly, the sizes measured in radio are not necessarily those of single merger components, but they should be interpreted primarily as separation between the two precoalescence starburst units (e.g., for all the three systems resolved in radio (ID 245158, 412250, 519651), their separation is larger than the size of single nuclei). Secondly, the early evolutionary phases are also characterized by lower dust obscurations, suggesting that the merger induced gas compaction (i.e., the increase of hydrogen column density in the center) has not yet completed.

This precoalescence subset identification provides an immediate physical interpretation for six galaxies of those that were simultaneously fitted with a double Gaussian in the 1D spectrum (Sect. 2.4), explaining this profile as coming from different merger components. However, we warn that these diagnostics are not identical and the connection between the line profiles in the 1D and in the 2D spectrum is not straightforward. Starbursts with multiple spatial emission lines do not necessarily display double Gaussians in the 1D spectrum, because this is subject to projection effects and depends on the distribution in wavelength of each spatial component. Indeed, the lines of one of the galaxies shown in Fig. 4 (ID 519651) were still fitted with a single Gaussian in the 1D.

Furthermore, our subset of six precoalescence starbursts identified from the 2D spectra is not necessarily complete, as many galaxies (e.g., ID 635862, 777034, 472775, 685067) have sky-subtracted 2D spectra with low S/N, not allowing to apply the visual criteria presented above in this Section. We would have required longer integration times or spatially resolved observations to build a complete sample of starbursts before the coalescence. Similarly, if the two merger nuclei are too close, it would be impossible to detect them even in the 2D spectra, and would need a significant improvement of spatial resolution to identify the pair.

4.2. Velocity enhancement and shocks toward the coalescence

4.2.1. BPT diagram and shocks

The second (bottom-left) panel of Fig. 3 shows that more obscured starbursts tend to have higher N2 relative to Hα, reaching [N II]/Hα ratios higher than 1, which are more typical of AGN and LINERs. Indeed, the classical BPT diagnostic diagram in Fig. 5, performed on nine galaxies with [O III] and Hβ available measurements, confirms that SBs with higher obscuration and line velocity width are found in the composite, AGN or LINER classification regions, according to empirical separation lines derived in the local Universe (Kauffmann et al. 2003; Kewley et al. 2001; Cid Fernandes et al. 2010; Veilleux & Osterbrock 1987).

Fig. 5. Top panel: BPT diagram of nine starbursts in our sample with optical spectra available (for one galaxy included in zCOSMOS, we did not detect both [O III]5007 and Hβ). While three sources lie in the SF excitation region, the remaining galaxies are not consistent with SF, and their spectra show a mixture of composite, AGN and LINER properties. The color coding indicates that galaxies with higher N2 which are closer to the AGN and LINER regions also have increasingly higher line velocity widths (σline). Bottom panel: same diagram as above, but here the galaxies are color coded according to their total dust attenuation AV, tot. More obscured starbursts preferentially display AGN or LINER properties.

Two variants of the BPT using the [SII]6717+6731/Hα or the [OIII]λ5007/[OII]λ3727+3729 ratios (S2BPT or O2BPT, respectively) are shown in Fig. B.1. We remind that, due to an enhanced ionization parameter and lower metallicity (at fixed mass) in the ISM at higher redshifts, the average star-forming galaxies population at z = 0.7 occupies a region in the BPT diagram which is shifted rightwards by ≲ + 0.1 dex compared to z = 0.1 (Faisst et al. 2018; Masters et al. 2016). However, there are currently no studies addressing how this will affect the separation lines among SB, AGN and LINERs. If we suppose that at z ∼ 0.7 the same shift applies also to these lines, galaxies at intermediate obscurations and line widths would still fall in the composite region with dominant LINER or AGN-like properties. Also, this would not affect our subsequent conclusions based on the comparison with shock models.

We noticed that the location of those galaxies shifted to the right compared to the purely SF region is consistent with the predictions of shock models, with varying shock contribution and velocity (compare with Figs. 10 and 2 of Rich et al. 2011, 2014, respectively). Additionally, Lutz et al. (1999) argue that LINER-like spectra in infrared selected galaxies are due to shocks, possibly related to galactic superwinds. The presence of increasing widespread shocks provides the most immediate interpretation for the spectra in our sample with enhanced [N II]/Hα, given that AGN emission would be highly suppressed (Paper I).

However, we cannot exclude some residual influence by an AGN. Hydrodynamical simulations performed by Roos et al. (2015) show that even in the case of high obscuration an AGN can ionize the gas very far from the nucleus, reaching kpc scales and the circum-galactic medium. Furthermore, the accreting black hole might not be in the center, but that sounds unplausible: the attenuations toward the center derived independently from the X-ray detected AGNs are consistent with those derived from the mixed model (see Paper I) and, even further, Rujopakarn et al. (2018) show that the AGN position correlates with that of active star forming regions. Finally, we also notice that two galaxies (which simultaneously have X-rays and mid-IR dusty torus detection) were fit with broad Hα components (line width of ∼1000 km s⁻¹). Such large velocity widths have been observed in both shock-dominated regions (possibly supernova driven, Ghavamian et al. 2017) and AGNs (Peterson 1997; Gaskell 2009; Netzer 2015). IFU data would be needed to disentangle shock or AGN emission, as we expect the latter to be much more concentrated in the central part of the system.

4.2.2. The dynamical masses of our sample

In order to better understand the dynamical status of our starbursts and how far they are from relaxation, we compared their total masses M_tot to the dynamical masses M_dyn estimated from the line velocity widths and radio sizes. For the latter, we used the formulation of Daddi et al. (2010) as :

$$\begin{array}{c}\hfill M_{\mathrm{dyn}}=1.3\times \frac{FWH{M}_{\mathrm{radio}}}{\times }{({\mathit{FWHM}}_{\mathrm{line},\mathrm{total}}/2)}^2{\mathrm{G}\sin }^2(i)\end{array}$$(4)

where FWHM_{line, total} is the one-dimensional total line velocity width (accounting for both rotation and dispersion), while sin²(i) is the correction for inclination that we take as the average value for randomly oriented galaxies (57°). In order to determine the total uncertainty on M_dyn, we considered an additional error on the inclination factor of 0.3 dex, as in Coogan et al. (2018). This represents the main contribution to the error (∼90% in median), since the line width and radio sizes are always well measured with high S/N.

Then we compared this quantity to the total mass content (baryonic + dark matter) of the systems, estimated as:

$$\begin{array}{c}\hfill M_{\mathrm{tot}}=M_{\ast }+M_{\mathrm{gas}}+M_{\mathrm{dark}\mathrm{matter}}\end{array}$$(5)

in which M_gas was determined, as described in Paper I, as M_gas = 8.05 + 0.81 × log(SFR_IR) (Sargent et al. 2014), valid for a starburst regime, and we assumed M_darkmatter = 10% ± 10% of M_*. Since this contribution is highly uncertain, it was set nearly uncontrained. However, this range is consistent with studies of high-z (>0.5) massive star-forming galaxies, which found a modest to negligible dark matter fraction inside the half-light radius (e.g., Daddi et al. 2010; Genzel et al. 2017). In any case, given the small contribution, its exact value does not affect the results of this paper. For the error determination, we considered the above uncertainty on M_darkmatter, a 0.1 dex error on M_* (Laigle et al. 2016), and 20% incertitude on the gas mass (even though its contribution is negligible given that M_gas ≃ 0.1 M_* on average for our sample).

The comparison between M_dyn and M_tot in Fig. 6 shows that, on average, our galaxies are not completely virialized: while nearly half of the sample is consistent within 2σ with the 1:1 relation, the remaining part is located above at higher M_dyn. The largest departures from virialization are observed for the precoalescence and less obscured systems, which are supposedly earlier stage mergers. On the contrary, the systems with better agreement may be fully coalesced starburst cores with higher A_V, tot.

Fig. 6. Top panel: comparison between the dynamical mass Mdyn and the total mass content Mtot (=Mgas + M* + Mdarkmatter) for our SBs sample, color coded by their total attenuation AV, tot. Bottom panel: diagram showing the square of the total FWHM velocity width as a function of Mtot/FWHMradio, using the same color coding based on AV, tot. On the y-axis, const = 1.3 × G/(4⟨sin2(i)⟩) groups the coefficients in Eq. (3) so as to facilitate comparison with the virialized case (1:1 relation, shown as a gray dotted line). The blue continuous line represents a linear fit to our sample, excluding galaxies with an upper limit on their radio size, while the blue shaded area shows the ±1σ limits of this best-fit relation. Both panels of the figure suggest that our galaxies may be approaching virialization, and more obscured starbursts are closer to the equilibrium.

The tight connection between velocity and gravitational potential is clarified in the bottom panel of Fig. 6, as the ${\text{FWHM}}_{\text{line}}^2$ and M_tot/FWHM_radio correlate at 5σ significance (with R = 0.74 and p-value = 0.0002). Also here, while precoalescence mergers have larger displacements from the 1:1 relation, they are confined in a region at lower velocity widths and shallower potential wells. This suggests that also other starbursts (ID 249989, 466112, 326384) in this region may be precoalescence mergers that we were not able to securely identify, due to their lower S/N 2D spectra, and indeed their optical morphology strengthens this suspicion. In the upper-right part of Fig. 6-bottom, separated from the previous sample, are clustered the more obscured starbursts, that is, supposedly coalesced mergers. We notice also that all X-ray detected AGNs are localized in this region of the diagram, indicating a possible link between evolutionary phase and AGN properties, that we will further investigate in the following Section.

Overall, the above results suggest a time-evolutionary scenario, in which more advanced, already coalesced mergers are close to virialization, and the increased central potential wells (due to the contribution of both merger components) are responsible for the enhancement of both the kinetic energy content and shocks toward later stages of the interaction. The tight relation between the line velocity width of single Gaussian components (a proxy for the velocity dispersion in the system) and shock production (traced by the N2 index) is further indicated by the color coding of the BPT diagram in Fig. 5-top, and by the correlation between FWHM_line and N II/Hα in Fig. 7, which has a significance higher than 5σ (R = 0.67, p-value = 3 × 10⁻⁴) and a dispersion of 0.26 dex.

Fig. 7. Correlation between the line velocity width of single Gaussian components (FWHMline) and the N2 index (both measured from our Magellan-FIRE spectra), indicating that the two quantities are tightly (scatter = 0.26 dex) physically related.

4.3. Lower line equivalent widths toward late merger stages

The equivalent widths (EW) of hydrogen recombination lines give a relatively dust-unbiased picture (assuming that stars and emission lines are equally extincted) of the contribution of the SFR to the stellar mass content, and they are sensitive to the luminosity weighted age of the stellar populations, so that they could provide useful information about the evolutionary stage of the merger. However, these EWs would only probe what is happening in the outer parts of the system, since the core is completely obscured in optical and near-IR.

We show in the last diagram in Fig. 3 and the upper part of Fig. 8 that, when the starbursts become more obscured, the EWs of Paβ, Hδ and Hβ decrease, indicating a gradual SFR decline in the outer skin of more obscured and compact starbursts. Additional correlations were found also independently between those EWs and the other quantities, such as the N2 index (Fig. 8-bottom). We additionally remark that the different Paschen and Balmer lines correlate each other (see Fig. B.3), for which reason our results, derived adopting the Paβ line as a reference (because it is the least attenuated), are also valid when considering the Hδ, Hβ and Hα lines.

Fig. 8. Top panels: comparison between the EW of Hδ, Hβ and Paβ lines to the dust attenuation parameter APaβ, IRX, defined as 2.5 × log10(1 + SFRIR/SFRPaβ, obs) (Paper I). We remark that the last panel is equivalent to the bottom-right plot in Fig. 3, even though a different scale has been used (the relation between APaβ, IRX and AV, tot is given in Eq. (1)). Bottom panels: correlations between the EW of Hδ, Hβ, Paβ lines and the N2 index (=log10([NII]/Hα)). The blue continuous lines are the best-fit linear relations, determined as explained in Sect. 3, while the blue shaded area show the ±1σ scatter of our data around the best-fit relations.

In our sample, we also found that the Balmer EWs, while having a large dynamical range, can reach very low values: in five galaxies (ID 303305, 685067, 777034, 862072 and 345018) we measured an EW(Hδ) <  −4 Å (i.e., in absorption), which are typically found in E+A dusty galaxies (Poggianti & Wu 2000). Low EW hydrogen recombination lines (in strong absorption) are clear signatures of the prevalence of A-type stars, indicating that a recent (<1 Gyr ago) massive star-formation episode has taken place during the past 10⁸−10⁹ yr, while the youngest stellar populations (mainly OB stars) are nearly all obscured by dust in the inner starburst core. Our dusty starburst systems should also not be confused with post-starburst (PSB) galaxies, which have similar absorption EWs (e.g., EW(Hδ) lower than ∼−5 Å as in the selection of Goto 2007 and Maltby et al. 2016), but are nearly or already quenched systems, owing much lower SFR levels and lower dust content compared to our sample (see Pawlik et al. 2018 for a full discussion of the different types of PSB galaxies). We caution that the quenched PSB selection from only the Balmer EW may be contaminated by real starburst systems.

Putting all together, the time-evolutionary scenario that we have suggested has the advantage of explaining in a simple way these new results. If we follow the merger evolution toward the coalescence, the outer starburst skin becomes increasingly dominated by A-type stars, recognizable through the deep absorption lines in the optical and near-IR, and which were formed at earlier times when the separation between the merging nuclei was larger. At the same time, the star-formation in the skin is being suppressed, possibly driven by supernova feedback.

4.4. Outliers

We found in Sect. 3 (Fig. 3) that four galaxies are outside the 1σ dispersion of the best-fit relations between the dust attenuation A_V, tot and, simultaneously, the [NII]/Hα ratio (N2), the line velocity width (FWHM_line) and the EW(Paβ). In particular, they have lower N2, FWHM_line and EW(Paβ) than expected from their A_V, tot, suggesting that, compared to other highly obscured galaxies, there is a minor impact from shocks or a dominant contribution of star-formation to the emission lines.

Within our SB sample, we recognize that these four outliers have the largest dust obscurations A_V, tot ≥ 18 mag, and are among the most compact, with radio FWHM sizes below 1 kpc. These extreme and peculiar features suggest they may represent the very end stages of the merger evolution, and that the correlations with A_V, tot may saturate toward these late phases. We also notice that the same objects are not systematically outliers when we consider their N2, FWHM_line and EW(Paβ) values, confirming the close physical connection among these quantities, as shown in Sects. 4.2 and 4.3.

4.5. The complete sequence of merger stages at intermediate redshift

Our observations and results, presented in previous sections, suggest we are starting to see an evolutionary sequence in high-redshift mergers. This can be traced through a variety of physical measurable quantities of our galaxies, including the total attenuation toward the center, the characteristic size of the starburst, the EW of hydrogen absorption lines, and finally the [N II]/Hα ratios and line velocity widths, which behave similarly. In Fig. 9 we schematize with a cartoon all the results that we have found so far, showing with a red continuous line the qualitative trend of the different physical quantities as a function of time. We divided the time axis into five merger evolutionary stages, which are arranged in relation to the two most crucial transformation events during the merger: the coalescence and the blow-out (or QSO) appearance.

Fig. 9. Schematic illustration of the time-evolutionary behavior of the physical parameters studied in the text: dust attenuation, characteristic size of the system, line velocity width (or, equivalently, the N2 index) and the EW of Balmer and Paschen lines. The time sequence is divided into five fundamental merger stages, with the QSO and passive spheroidal system representing the final stages according to the classical merger paradigm (Sanders & Mirabel 1996; Hopkins et al. 2008a,b). Solid lines are qualitative trends during the SB phase inferred from our results, while dashed lines are predictions for the future evolution of the four parameters shown on the y-axis (line width and N2 index behave similarly). In the upper part of the figure, we show a qualitative merger timescale following Fig. 1 of Hopkins et al. (2008b), assuming for the merger a total starburst duration of 200 Myr. For each phase, we show in the bottom part the ACS-F814W cutout of a representative case. The first three images are SB galaxies from our sample: ID 223715, ID 777034 and ID 472775. They were chosen as having increasing dust attenuations and radio compactness, suggestive of more advanced merger phases: the first was identified as a precoalescence merger in Sect. 4.1, while the latter is unresolved in radio and is highly obscured (AV, tot = 18 mag). The last two cutouts show a quasar at z = 0.73 and an ETG at z = 0.66, selected in COSMOS field from the catalogs of Prescott et al. (2006) and Tasca et al. (2009), respectively.

We notice that the first phase may not necessarily represent the beginning of the interaction, that is, when the two galaxies approach for the first time. Even though the whole merger episode may last 1–1.5 Gyr in total, from the first encounter to the formation of a passive spheroidal system, the starburst activity is typically shorter, ranging 200–300 Myr (Di Matteo et al. 2008) and may be triggered intermittently at various stages of the evolution. Furthermore, whether or not a strong burst is already activated at the first approach depends on many factors, including the impact geometry, the morphology, the stellar mass ratio and the gas content of the colliding galaxies (Di Matteo et al. 2008).

Besides the observable starburst phases studied in this work, can we also make some predictions on the future evolution of these systems? In general, it is very hard to demonstrate visually a connection between mergers and their descendants. Indeed, not all merger-induced starbursts exhibit morphological disturbances (Lotz et al. 2008), and when merger residual signatures are present, they fade rapidly, becoming almost invisible beyond the local Universe even in the deepest optical images (e.g., Hibbard & Vacca 1997). We can in principle rely on hydrodynamical simulations, which allow to trace the full time-sequence of mergers, even though they also present limitations due to the many assumptions, initial conditions and physical complexity involved in such events.

In the classical theoretical merger paradigm, the infalling gas triggers obscured AGN accretion (Bennert et al. 2008a), whose peak of activity typically occurs ∼250 Myr after the onset of the starburst (Wild et al. 2010), and ∼100 Myr after the peak of SFR (Davies et al. 2007; Hopkins 2012). It is during these later starburst phases that the AGN feedback can blow out with strong feedback winds the surrounding dust and gas cocoon, eventually revealing itself as a bright QSO (Hopkins et al. 2008b). This phase is generally very short, lasting for ≲100 Myr (Hopkins & Elvis 2010), and has been claimed since a long time: Lípari et al. (2003) suggested that QSOs could be indeed young IR active galaxies at the end phase of a strong starburst.

Since the QSO dominates the luminosity of the system at all wavelengths, it would be extremely hard to analyze the physical properties of the host galaxies during this phase. Indeed, Zakamska et al. (2016) show that even in radio-quiet QSOs both the infrared and the radio emission are dominated by the quasar activity, not by the host galaxy. An alternative possibility is to look far from the central bright source. Recent works are revealing Ly-α nebulae surrounding high-redshift quasars, with extension that can reach tens of kpc (≲50 kpc) from the center (Arrigoni Battaia et al. 2019). On the other hand, one may focus on local samples, increasing simultaneously the images resolution. For example, Lípari et al. (2003) and Bennert et al. (2008b) discovered with HST the presence of outflows, arcs, bubbles and tidal tales in optical band in a sample of local QSOs, possibly formed through strong galactic winds or merger processes. Again in nearby (z <  0.3) QSOs, near-IR H band adaptive optics observations (Guyon et al. 2006) revealed that ∼30% of their hosts show signs of disturbances, and the most luminous QSOs are harbored exclusively in ellipticals or in mergers (which may become ellipticals soon). Furthermore, while the SFRs of the hosts are similar to those of normal star-forming galaxies, their mid- and far-IR colors resemble those of warm ULIRGs, strengthening a connection between these two objects.

In the following two Sects. 4.6 and 4.7, we discuss separately the two ending stages of the merger sequence, and investigate how our work can provide some clues to understand what are the physical properties of the systems into which our starbursts will evolve. In the cartoon of Fig. 9, the predicted evolution for all the quantities studied in this paper (see Sect. 3) is shown with a dashed line. These qualitative trends are motivated mainly from simulations, and are not confirmed observationally.

4.6. Mass-size relation and comparison with higher and lower-z starbursts

The merger-induced starbursts are supposed to end up in a passive system, but we do not know the exact physical properties (e.g., size, stellar mass, morphology) of these merger remnants. Sub-millimeter galaxies (SMGs) at high redshift (>2), which are commonly viewed as higher luminous counterparts of lower redshift ULIRGs, have been suggested to be direct progenitors of massive ETGs (Tacconi et al. 2008; Toft et al. 2014). We can investigate this connection by comparing in Fig. 10 the stellar masses and the characteristic sizes of our starbursts with those of disk galaxies and spheroids at z ∼ 0.7 (van der Wel et al. 2014). To be conservative, we are adopting here the M_*-size relations for circularized radii. If we consider instead the non-circularized cases, the same relations would slightly shift upwards by ∼0.1 dex. In addition, we remark that we are comparing our radio (starburst) extensions to optical rest-frame sizes tracing the stellar mass distribution of disks and elliptical galaxies. Indeed, we implicitly assume that, after the gas in our starburst cores is converted into stars, the extensions of these cores will represent also the stellar component sizes of their passive remnants. On the other hand, they may still represent the dense star-forming gas components of post-starburst systems if some residual is left after the merger, as they may remain compact for at least 1 Gyr (Davies et al. 2019).

Fig. 10. Diagram showing the radio size vs. stellar mass for our sample. We compare our results to the stellar mass–stellar size relation of LTGs (cyan line with 1σ dispersion) and ETGs (blue line) at z = 0.7 (van der Wel et al. 2014), and with the bulge properties of low redshift (z ∼ 0.1) spiral galaxies from Graham & Worley (2008; gray squares). The eleven points shown here for the bulges represent the median (±1σ) of their distributions of stellar masses and stellar sizes (in K band) as a function of galaxy type, from S0 to Sm spirals. For three galaxies in our sample fit with a double Gaussian, we also represent the “deblended” radio sizes of each single component with violet crosses, connecting them with a dashed violet line. In these cases, we assigned to each component half of the total stellar mass of the system, even though a precise estimation requires a separate fit on deblended photometric data.

In the diagram of Fig. 10, six SBs are consistent with the late-type galaxy (LTG) relation at z ∼ 0.7. However all of them are precoalescence SBs and, as we have seen before, they should not be considered disk galaxies as their size is primarily reflecting the separation between the merging components. For two of the three galaxies resolved in radio, the single values return below on the early-type galaxy (ETG) relation. The characteristic sizes of this subset (FWHM_size ranging 3–15 kpc in diameter, with median FWHM_size of 8 kpc) are similar to those typical observed in SMGs (Casey et al. 2011; Tacconi et al. 2008; Biggs & Ivison 2008), which suggests that SMGs at high-redshift (or at least a fraction of them) could be indeed intermediate-phase mergers composed of unresolved double nuclei, as argued by Iono et al. (2009) and Arribas et al. (2012).

In the bottom part instead, we can immediately notice that a major fraction of our sample (13 galaxies, i.e., 52% of the total sample) is not consistent with the ETG relation (taking 1σ dispersion), and is located well below it by ∼0.5 dex, with an average size of ≤1.2 kpc, indicating that they are much more compact than their stellar envelopes and than typical ellipticals at z ∼ 0.7. We underline that such difference would be even higher if we compare this subset to the M_*-size relation at redshifts lower than 0.7, as the ETG sizes at z = 0.25 are a factor of 1.5 higher than those at z = 0.75, at our median stellar mass (van der Wel et al. 2014). This sample of very compact starbursts has typical extensions that are similar to those of dense star-forming regions in local ULIRGs (Genzel et al. 1998; Piqueras López et al. 2016), including Arp 220 (Sakamoto et al. 2017) and M82 (Barker et al. 2008), suggesting they are driven by the same merger mechanisms (as also argued in Paper I).

If we take for each galaxy its distance from the LTG relation (dist_LTG = log₁₀(FWHM_size/FWHM_LTG)), we can also use this quantity in place of the radius to trace the same sequence found in Sect. 3, taking into account the mild dependence on stellar mass. As the merger proceeds, the system moves from the LTG to the ETG relation and then even below at significantly smaller sizes (by ∼0.5 dex at least), meaning that the compact starburst cores that form at the coalescence cannot produce directly the ellipticals seen at redshift 0.7 and below.

The sizes of our starbursts instead resemble those of typical bulges in lower redshift spirals and lenticular galaxies (Graham & Worley 2008; Laurikainen et al. 2010), indicating a possible evolutionary link between mergers and bulges, as suggested by other works (e.g., Sanders & Mirabel 1996; Lilly et al. 1999; Eliche-Moral et al. 2006; Querejeta et al. 2015). This idea is consistent with the typical observed gas fractions of our starbursts (derived as M_gas/(M_* + M_gas), with M_gas calculated in Sect. 4.2.2), which range between 0.02 and 0.25 (∼0.1 in median). Assuming that all the remaining gas is consumed before the passivization and that the same amount of gas has been already converted into stars (which depends on the merger phase and dynamics), it means that the current starburst cores can produce approximately 20%, and up to 50%, of the final stellar mass of the galaxies. Higher resolution radio images targeting specific emission lines can further constrain the kinematic properties of the starbursting cores, by looking for rotation, or their luminosity profile, for example measuring their Sersic index. How this old stellar component is affected by the merger can depend on many conditions difficult to model in detail, including the geometry of the interaction, the gas content and the mass ratio of the colliding galaxies.

4.7. QSOs in formation at z ∼ 0.7?

In the starburst selection phase, we discarded several quasars because of the impossibility to study the properties of their host galaxies (dust attenuation, SFR, stellar mass), as discussed in Sect. 4.5. In order to overcome these limitations, several authors have extensively studied also the transitional moments (just preceding the final blow-out) in the dress of type-I and warm ULIRGs (Kawakatu et al. 2006; Sanders et al. 1988). Similarly, we can have some clues of the inner black-hole activity just before the hypothetical QSO by looking at the AGN diagnostics in our starburst sample.

As mentioned in Sect. 2.8, we detected the mid-IR dusty torus AGN emission in 12 galaxies, and simultaneous X-ray emission in six of them. Notably, the latter are the only ones (among mid-IR AGNs) whose host galaxies lie below the ETG relation in the mass-size diagram (see Fig. 10), at systematically smaller sizes than ellipticals at z ≤ 0.7. This suggests that during early merger stages the AGNs are predominantly obscured, while they start to appear in X-rays toward intermediate stages (i.e., when the host starbursts are more compact and obscured), possibly driven by rapid AGN feedback clearing the gas and dust around the black hole.

This can be seen better in Fig. 11-left, where all the host galaxies of X-ray detected AGNs are located at larger distances from the mass-size relation of LTGs compared to X-ray undetected AGNs. Moreover, they have X-ray luminosities at least 1 order of magnitude higher (∼ + 1.5 on average) than what inferred from their SFR. For the X-undetected galaxies instead, the upper limit on L_X = 10^41.7 erg s⁻¹, determined by average-stacking their fluxes in the 2−10 keV band at the median z of the sample, is consistent with emission produced by star-formation only, suggesting that in this band the AGN is completely obscured.

Fig. 11. Left panel: comparison between the observed X-ray luminosity (LX, obs) and the distance from the mass-size relation of LTGs at z ∼ 0.7, for our six starbursts detected in X-rays. Upper limit on LX, obs for six mid-IR AGNs undetected in X is shown with a black circle, where the horizontal segment represents the range of distLTG spanned by this subset. The intrinsic X-ray luminosity due to star-formation is highlighted with a gray line for the median SFR of the sample (±0.4 dex scatter from Mineo et al. 2014), and may dominate the total X-ray observed emission for the X-undetected starbursts. Center panel: X-ray attenuation LX, obs/LX, int as a function of the infrared-based attenuation AV, tot (in a mixed model geometry and toward the center) for our sample of mid-IR detected AGNs. We assumed here a bolometric correction factor LX, intr, AGN = 0.04 × LBOL, AGN (Vasudevan & Fabian 2007). Stacks on the whole sample and on the X-undetected subset are displayed with hatching circles, while the violet shaded regions indicate the area of no obscuration, which incorporates a factor of two uncertainty in the conversion between intrinsic X-ray and bolometric AGN luminosity. Right panel: same diagram as before, but assuming an LBOL, AGN- dependent bolometric correction (Lusso et al. 2012), as explained in the text.

In order to assess the level of obscuration, we computed the ratio between the observed and the intrinsic X-ray luminosity L_X, obs/L_X, intr, comparing this quantity to the total dust attenuation A_V, tot inferred in a mixed model from Paβ and the bolometric IR luminosity (Paper I). L_X, intr comprises the contribution from both star-formation (as explained in Sect. 2.8) and from the AGN, assuming L_{X, intr, AGN} = 0.04 × L_BOL, AGN (Vasudevan & Fabian 2007) and a bolometric AGN luminosity L_BOL, AGN = 1.5 × L_AGN, IR (Elvis et al. 1994; Fig. 11-center). Alternatively, we considered the bolometric correction of Lusso et al. (2012) for type-2 AGNs, which depends on L_BOL, AGN itself through the following, nonlinear equation:

$$\begin{array}{c}\hfill {log}_{10}\left(\frac{L_{\mathrm{BOL}}}{L_{2-10\mathrm{keV}}}\right)=0.217x-0.022x^2-0.027x^3+1.289,\end{array}$$(6)

where x = log₁₀(L_BOL)−12 and the scatter of the relation is 0.26 dex. However, the results derived with this second assumption (Fig. 11-left) do not change significantly compared to the first case.

The figures presented above confirm that the total attenuation inferred from L_X may be the discriminating parameter between X-ray detected and undetected mid-IR AGNs: while the former are relatively unobscured, the X-ray emission from the second is suppressed at least by a factor of 30. Interestingly, this transition does not seem to be related to an increased bolometric luminosity, since no correlations are observed with this quantity. We also remind that for this analysis we are following the standard procedure, which does not take into account shock contribution to the X-ray luminosity. A possible non negligible shock emission at these high energies (which in any case is difficult to model) would result in an underestimation of the true effective X-ray attenuation toward the AGN.

The previous results suggest that the X-ray attenuation decreases as the starburst becomes more dust-obscured (probed by A_V, tot) during the last merger phases. In a standard framework (i.e., if we exclude dominant contribution from shocks to the X-ray luminosity), this apparent contradiction can be reconciled by considering the different timescales of our diagnostics. On the one hand, Paβ (used to calculate A_V, tot), yields a luminosity (or, equivalently, a SFR, by applying the Kennicutt et al. (1994) conversion) that is averaged over a timescale of 20−30 Myr. Conversely, the AGN luminosity that we measure in X-rays gives an istantaneous information of the AGN activity. As a consequence, with the X-ray analysis we are able to probe the current dust attenuation level, while A_V, tot traces the obscuration in the recent past (≲30 Myr). According to this speculation, the X-ray luminosities measured for a subset of six late stage mergers indicate that the AGN-induced blow-out may have already started since a few Myr ago, clearing the surrounding gas and dust content, and that we might be very close to the final QSO phase.

Furthermore, in Fig. 12 we display the AGN accretion efficiencies of our galaxies as a function of their distance from the LTGs mass-size relation. The efficiencies were estimated by comparing the observed L_BOL, AGN/M_* ratios to the maximum value allowed by Eddington (L_BOL, AGN/M_*|_EDD ≃ 1.5), from which we derived the so-called Eddington ratio (L/L_EDD)|_AGN. We assumed the typical correlation for AGNs between the stellar mass M_* and black hole mass M_BH of Reines & Volonteri (2015) and a spherically symmetric accreting BH, yielding log(L_EDD/M_*)≃0.9685 + 0.05log₁₀(M_*). The M_* in the second term can be approximated with the median value of the sample M_*, median, leaving a small secondary dependence on stellar mass which, for our mass ranges (10¹⁰ − 10¹¹ M_⊙) produces variations of <5%. This variation, added to the uncertainty on the relation between M_* and M_BH reported by Reines & Volonteri (2015) (cf. their Eqs. (4) and (5)), is highlighted in Fig. 12 with a blue shaded area around the Eddington limit (blue line) calculated above.

Fig. 12. LBOL, AGN/M* vs. distance from the mass-size relation of LTGs (distLTG) for our SBs sample. The Eddington limit is shown with a blue horizontal line, while the shaded area takes into account the spread of M* among our sample and the uncertainty of the relation between stellar mass and BH mass by Reines & Volonteri (2015). The Eddington ratio is shown on the right y-axis.

From this analysis, we found that two AGNs have an Eddington ratio higher than 1 (1.2 < (L/L_EDD)|_AGN <  1.85), while additional three AGNs are radiating between 57% and 79% of their maximum luminosity. However, all these five AGNs are still consistent within 2σ errors with (L/L_EDD)|_AGN = 1 if we also consider the uncertainty on the conversion factor from the ratio L_BOL, AGN/M_*, as discussed before. We remark that additional uncertainties on the M_* − M_BH relation, which depends on the assumptions about the BH accretion geometry, are not taken into account here. The remaining seven IR-detected AGNs have instead lower Eddington ratios between 0.35 and 0.08.

The galaxies which are undetected in X, radio and mid-IR may contain low-active AGNs, even though current upper limits on the Eddington ratio are not so stringent and do not allow to discriminate them from the detected subset. The intrinsic variability of AGN accretion may thus explain why we are currently missing many of these sources in our sample, and that only deeper X-ray observations can potentially reveal. The duty cycles above 30% and 1% L_EDD seem to be at least ∼25% and ∼50%, respectively.

5. Summary and conclusions

Using our unique sample of 25 starburst galaxies (typically seven times above the star-forming main sequence) at z = 0.5−0.9 with near-IR rest frame spectroscopy of Paschen lines, we found in Paper I that they span a large range of attenuations toward the core centers from A_V = 2 to A_V = 30, forming a sequence which is consistent with a mixed model geometry of dust and stars. In this paper we have investigated the nature of this attenuation sequence, comparing A_V with other physical properties, such as the radio size (which traces the extension of the starburst), the emission lines velocity widths and [N II]/Hα ratios (which depend on the internal motions of the gas and the contribution from shocks), and finally the EW of hydrogen absorption lines, which is sensitive to the luminosity-weighted age of the stellar populations surrounding the optically obscured core.

We summarize the main results of this paper as follows:

1. We found that the physical quantities introduced above, namely the radio sizes (FWHM_radio), the line velocity widths (FWHM_line), the [N2 ]/Hα ratios (N2) and the equivalent widths of Paschen and Balmer lines (EW_{Balmer, Paschen}), all correlate each other (Fig. 3), defining a one-parameter sequence of z ∼ 0.7 starburst galaxies.

2. These correlations can be interpreted as a time-evolutionary sequence of merger stages. As the merger evolves, the starburst becomes more compact and dust obscured, while the deep potential wells created by merging nuclei produce, according to the virial theorem, an increase of the kinetic energy and shocks in the system. At the same time, intermediate aged A-type stars in the outer starburst core regions are primarily responsible for the stronger optical+near-IR absorption lines in later phases.

3. Four galaxies are outliers simultaneously in three of the ten main correlations, which involve A_V, tot and, respectively, N2, FWHM_line and EW(Paβ). Having the largest dust attenuations and among the smallest radio sizes in our sample, these outliers may represent the very end phases of the merger evolution, where the above three relations may reach a saturation level.

4. Using sky-subtracted 2D spectra, we identified a subset of seven precoalescence mergers by the presence of spatially separated or kinematically detached (i.e., rotation-driven tilted lines with different inclination angles) Hα components, representing earlier, less obscured phases of the interaction. The radio sizes measured for these systems are likely tracing the separation between the merging nuclei rather than the dimensions of single cores. However, our sample may contain additional double nuclei which we are not currently able to resolve.

5. Half of our sample comprises extremely compact starbursts, with average half-light radii of 600 pc (six galaxies have only upper limits), similar to the sizes of starbursting cores observed in local ULIRGs. These sizes are also ∼0.5 dex smaller than ETGs at redshift ∼0.7 and below, indicating that our merger-driven starbursts cannot be direct progenitors of the population of massive ellipticals formed in the last ∼7 Gyr. On the contrary, they are more consistent with typical sizes and masses of bulge structures in spiral galaxies (Graham & Worley 2008), suggesting a possible evolutionary connection between our starburst cores and bulges.

6. In our sample, we detected at >3σ the mid-IR dusty torus AGN emission in 12 starbursts, with Eddington ratios ranging from 1.9 to <0.08. Among them, only six galaxies are simultaneously detected (at 3σ) in X-rays. Intriguingly, the latter have the largest departures from the mass-size relation of LTGs (at z ≤ 0.7), suggesting that AGNs start to appear in X-rays during the latest (compact) merger phases, as the blow-out of surrounding dust and gas may precede a possible final QSO.

Overall, the relations among the above physical parameters converge toward a time-evolutionary sequence of merger stages, which represents an observational evidence (translated at higher redshift) of the theoretical merger-induced starbursts framework of Hopkins et al. (2008a,b), Di Matteo et al. (2005), and the evolutionary sequence postulated by Toomre & Toomre (1972). The future advent of JWST will allow to test this scenario up to very high redshift in different galactic physical conditions compared to the epochs studied here.

Acknowledgments

We thank the anonymous referee for useful suggestions that improved the quality of this manuscript, G.Rudie for assistence with Magellan observations, Nicolás Ignacio Godoy for data reduction, and Daniela Calzetti for discussions. M.O. acknowledges support from JSPS KAKENHI Grant Number JP17K14257. N.A. acknowledges support from the Brain Pool Program, funded by the Ministry of Science and ICT through the Korean National Research Foundation (2018H1D3A2000902). M.B. acknowledges FONDECYT regular grant 1170618. R.C. acknowledges financial support from CONICYT Doctorado Nacional No. 21161487 and CONICYT PIA. ACT172033. A.C. acknowledges RadioNet conference funding. This research has made use of the zCosmos database, operated at CeSAM/LAM, Marseille, France.

References

Appendix A: Additional table

Appendix B: Additional plots

We display in Fig. B.1 two different versions of the BPT diagram, involving the [S II]6717+6731/Hα line ratio (top panel), and the [O III]5007/[O II]3727+3729 ratio (bottom panel), color coded according to the total velocity width of the lines (σ_line).

Fig. B.1. Diagram of [O III]/Hβ compared to [S II]6717+6731/Hα (top) and [O III]5007/[O II]3727+3729 vs. [N II]6583/Hα (bottom) for 8 galaxies of our sample with optical spectra available and contemporary detection of [O III] and one between [S II]+Hβ or [O II]. The two diagrams are also called the S2 and the O2 BPT diagrams, respectively. As in the classical BPT diagram in Fig. 5, they show that our starbursts have different line excitation properties, and those with higher line velocity widths are generally shifted toward the AGN and LINER regions, according to empirical separation lines derived in the local (z ≲ 0.3) Universe.

In Fig. B.2 we show the SED fitting of all the 25 starbursts in our sample, performed as explained in Sect. 2.1. The majority of our galaxies were fit with a GN20 template from Magdis et al. (2012). However, for six objects (ID 245158, 223715, 470239, 493881, 635862 and 862072), a MS template gave a better χ² than the GN20 SED. This is due to the fact that this particular template is not universal for all starbursts, and the true SB SED can have multiple dust temperatures.

Fig. B.2. SED fitting of the 25 starburst galaxies studied in this paper, presented in numerical order. The SED fitting procedure is described in Sect. 2.1. Blue and red curves show the stellar component (Bruzual & Charlot 2003) and the AGN torus emission (Mullaney et al. 2011), while the dust continuum emission is shown in green (Magdis et al. 2012). Downward arrows show the 2σ upper limit photometry at given wavelength. In each panel, zspec is the spectroscopic redshift from our Magellan spectra, zphot is the photometric redshift from Laigle et al. (2016), z is the output redshift of the best-fit SED (which in our cases was fixed to the spectroscopic value), and SFR is the best-fit infrared star-formation rate. The vertical dotted lines indicates the optimized range for fitting the stellar SED only (in order to avoid contamination from AGN torus and dust emission). The rightmost line also indicates the starting wavelength for fitting the AGN + dust SED in our paper. We remind the reader that the radio photometry (at 1.4 and 3 GHz) was not used in the fit, as explained in Sect. 2.1.

Fig. B.2. continued.

Fig. B.2. continued.

Finally, we display in Fig. B.3 the correlations among Hδ, Hβ, Hα and Paβ equivalent widths.

Fig. B.3. Comparison between the EWs of Hδ, Hβ, Hα and Paβ lines. The blue lines represent the best-fit linear relation, while the blue shaded areas indicate the ±1σ dispersion of the galaxies around the best-fit line. We display in top-left corner of each panel the Spearman correlation coefficient “R” and corresponding “p”-value, the residual χ2 of the correlation (“${\chi }_{\text{red}}^2$”) and the 1σ scatter around the best-fit line, whose equation (comprising 1σ errors on the best-fit coefficients) is shown in the bottom-right corner. In the four diagrams, all the points with available measurements (excluding upper limits) were used in the fit.

Appendix C: Hα and Paβ emission line fits

In this section, we show in Figs. C.1, the Hα+[N II]λ6583 Å and Paβ emission lines (if detected) coming from the Magellan-FIRE spectra of our starbursts sample, fit with single or double Gaussian profiles by using the tool MPFIT (see Sect. 2.4 for the procedure). We show the best-fit Gaussians with red lines, superimposed on the fully reduced and calibrated spectra (black line). The noise spectrum in the same spectral range of the object is drawn in the lower panel for each cutout (green line), while the fit residuals (noise normalized, as described in Calabrò et al. 2018) are displayed in the bottom panel with a black line.

For 12 galaxies in our sample, we fit their emission lines with a double Gaussian function (see Sect. 2.4). Among them, six are also precoalescence SBs (ID 245158, 493881, 470239, 223715, 545185, 668738) (see Sect. 4.1), and we interpret the two 1D Gaussians as coming from different merger components.

For the starbursts ID 862072 and 249989, the shape of the Hα line in the 2D spectra is consistent with global rotation, as they show a spatially extended Hα line with a single, uniform inclination. Additionally, their 1D line profiles are nearly symmetric, that is, the two Gaussian components have the same flux within the uncertainties. Even so, we cannot definitely exclude that these signatures are also produced by identical merging pairs rotating around the common barycenter with opposite relative velocities.

In four cases (ID 635862, 777034, 472775, 685067) we were not able to apply the two visual criteria presented in Sect. 4.1, because of the poor S/N of their sky-subtracted 2D spectra. However, we noticed that three of them have very asymmetric 1D emission line profiles, with fluxes of the two Gaussian components differing by >50%. Since they are less likely to be produced by a single rotating disk, we suggest they might be due to different merger components.

Fig. C.1. For each galaxy, in two contiguous panels (top and bottom): Hα+[N II]λ6583 Å and Paβ emission lines from the Magellan-FIRE spectrum (Paγ is shown for the galaxy ID 245158, since this line is used to infer the attenuation). Each panel is further divided into three sections: in the upper part is shown with a black line the 1D portion of the spectrum (fully calibrated) close to Hα or Paβ line, along with the best-fit Gaussian superimposed (red line). In the second part below we display the noise spectrum for the same spectral range of the object (green line), while the fit residuals (noise normalized, as described in Calabrò et al. 2018) are shown in the bottom section with a black line. The galaxies are presented in numerical order for double detections, while those with only Hα+[N II] detection are shown at the end.

Fig. C.1. continued.

Fig. C.1. continued.

Fig. C.1. continued.

All Tables

All Figures

	Fig. 1. For each of the galaxies observed in our second observing run at Magellan from left to right: HST-ACS F814W images (FWHMres = 0.095″), H-band UltraVISTA cutouts (same f.o.v. and FWHMres ∼ 0.75″) and 3 GHz radio images from VLA-COSMOS 3 GHz large project (Smolčić et al. 2017; FWHM ∼ 0.75″). Cutout images for the sample of the first observing run were presented in Paper I.
In the text

Fig. 2. GALFIT fitting of the radio VLA image (3 GHz) of the galaxy ID 685067 (z = 0.37) with the VLA synthesized beam in the upper row and with a Gaussian profile (convolved with the beam) in the bottom row. In horizontal sequence are shown, from left to right, the original image, the fitted model and the residual (original-model). For this galaxy, we derive an angular size of 0.20 ± 0.04″ (the pixel scale is 0.2″ pixel−1), which corresponds to a physical size of 1.06 ± 0.19 kpc. We notice that here GALFIT converges when fixing the position angle (PA) and axis ratio (q) parameters, to 0 and 1, respectively (Table A.1). This example illustrates the possibility to reliably measure the radio sizes of our objects even when they are smaller than the FWHM resolution of VLA (0.75″). In this case, the difference between the two models is recognized by looking at the residual images (i.e., original-model).

In the text

Fig. 3. Left column: correlations of AV, tot with the radio size (top), and with the N2 index (bottom). We show with filled blue diamonds all the Magellan SBs that were used to derive the best-fit linear relation (blue continuous line) and the ±1σ dispersion (blue shaded area), while empty diamonds represent galaxies excluded from those calculations. The latter comprise the four outliers discussed in the text (ID 303305, 345018, 472775, 545185) and all the upper limits in the AV, tot vs. FWHMradio plot. In the corners, we show in black the equation of the linear fit (which includes 1σ error of the two best-fit parameters), and in gray the Spearman correlation coefficient (R) with the corresponding p-value (p), the reduced chi-square of the fit (${\chi }_{\text{red}}^2$), and the 1σ scatter of our SBs around the best-fit line, all of which do not take into account upper limits and the four outliers mentioned above. For comparison, the linear fit and 1σ dispersion including the four outlier galaxies are highlighted with a gray continuous dashed line and two dotted lines of the same color. Right column: correlations of AV, tot with the line velocity width (top), and with the EW(Paβ) (bottom). In the last diagram, four galaxies without EW(Paβ) measurements are not considered.

In the text

Fig. 4. Close-up view of the Hα emission line profiles for the galaxies satisfying one of the two precoalescence criteria defined in Sect. 4.1. In each panel, as a consequence of the sky-subtraction procedure applied to the 2D spectra, the lines appear twice in different slit positions, in the first with a positive flux (in yellow) and in the second with a negative flux (in black). For each cutout, we show with two white arrows the slit position (s) and the dispersion direction (λ), which are slightly rotated due to the curved Magellan spectra. In order to clarify the classification criteria adopted in this work for finding precoalescence SBs, we highlight: (1) with continuous lines the different tilting angles of Hα line profiles (first three SBs); (2) with dotted ellipses the spatially separated Hα lines coming from different merger components (last four SBs).

In the text

Fig. 5. Top panel: BPT diagram of nine starbursts in our sample with optical spectra available (for one galaxy included in zCOSMOS, we did not detect both [O III]5007 and Hβ). While three sources lie in the SF excitation region, the remaining galaxies are not consistent with SF, and their spectra show a mixture of composite, AGN and LINER properties. The color coding indicates that galaxies with higher N2 which are closer to the AGN and LINER regions also have increasingly higher line velocity widths (σline). Bottom panel: same diagram as above, but here the galaxies are color coded according to their total dust attenuation AV, tot. More obscured starbursts preferentially display AGN or LINER properties.

In the text

Fig. 6. Top panel: comparison between the dynamical mass Mdyn and the total mass content Mtot (=Mgas + M* + Mdarkmatter) for our SBs sample, color coded by their total attenuation AV, tot. Bottom panel: diagram showing the square of the total FWHM velocity width as a function of Mtot/FWHMradio, using the same color coding based on AV, tot. On the y-axis, const = 1.3 × G/(4⟨sin2(i)⟩) groups the coefficients in Eq. (3) so as to facilitate comparison with the virialized case (1:1 relation, shown as a gray dotted line). The blue continuous line represents a linear fit to our sample, excluding galaxies with an upper limit on their radio size, while the blue shaded area shows the ±1σ limits of this best-fit relation. Both panels of the figure suggest that our galaxies may be approaching virialization, and more obscured starbursts are closer to the equilibrium.

In the text

	Fig. 7. Correlation between the line velocity width of single Gaussian components (FWHMline) and the N2 index (both measured from our Magellan-FIRE spectra), indicating that the two quantities are tightly (scatter = 0.26 dex) physically related.
In the text

Fig. 8. Top panels: comparison between the EW of Hδ, Hβ and Paβ lines to the dust attenuation parameter APaβ, IRX, defined as 2.5 × log10(1 + SFRIR/SFRPaβ, obs) (Paper I). We remark that the last panel is equivalent to the bottom-right plot in Fig. 3, even though a different scale has been used (the relation between APaβ, IRX and AV, tot is given in Eq. (1)). Bottom panels: correlations between the EW of Hδ, Hβ, Paβ lines and the N2 index (=log10([NII]/Hα)). The blue continuous lines are the best-fit linear relations, determined as explained in Sect. 3, while the blue shaded area show the ±1σ scatter of our data around the best-fit relations.

In the text

Fig. 9. Schematic illustration of the time-evolutionary behavior of the physical parameters studied in the text: dust attenuation, characteristic size of the system, line velocity width (or, equivalently, the N2 index) and the EW of Balmer and Paschen lines. The time sequence is divided into five fundamental merger stages, with the QSO and passive spheroidal system representing the final stages according to the classical merger paradigm (Sanders & Mirabel 1996; Hopkins et al. 2008a,b). Solid lines are qualitative trends during the SB phase inferred from our results, while dashed lines are predictions for the future evolution of the four parameters shown on the y-axis (line width and N2 index behave similarly). In the upper part of the figure, we show a qualitative merger timescale following Fig. 1 of Hopkins et al. (2008b), assuming for the merger a total starburst duration of 200 Myr. For each phase, we show in the bottom part the ACS-F814W cutout of a representative case. The first three images are SB galaxies from our sample: ID 223715, ID 777034 and ID 472775. They were chosen as having increasing dust attenuations and radio compactness, suggestive of more advanced merger phases: the first was identified as a precoalescence merger in Sect. 4.1, while the latter is unresolved in radio and is highly obscured (AV, tot = 18 mag). The last two cutouts show a quasar at z = 0.73 and an ETG at z = 0.66, selected in COSMOS field from the catalogs of Prescott et al. (2006) and Tasca et al. (2009), respectively.

In the text

Fig. 10. Diagram showing the radio size vs. stellar mass for our sample. We compare our results to the stellar mass–stellar size relation of LTGs (cyan line with 1σ dispersion) and ETGs (blue line) at z = 0.7 (van der Wel et al. 2014), and with the bulge properties of low redshift (z ∼ 0.1) spiral galaxies from Graham & Worley (2008; gray squares). The eleven points shown here for the bulges represent the median (±1σ) of their distributions of stellar masses and stellar sizes (in K band) as a function of galaxy type, from S0 to Sm spirals. For three galaxies in our sample fit with a double Gaussian, we also represent the “deblended” radio sizes of each single component with violet crosses, connecting them with a dashed violet line. In these cases, we assigned to each component half of the total stellar mass of the system, even though a precise estimation requires a separate fit on deblended photometric data.

In the text

Fig. 11. Left panel: comparison between the observed X-ray luminosity (LX, obs) and the distance from the mass-size relation of LTGs at z ∼ 0.7, for our six starbursts detected in X-rays. Upper limit on LX, obs for six mid-IR AGNs undetected in X is shown with a black circle, where the horizontal segment represents the range of distLTG spanned by this subset. The intrinsic X-ray luminosity due to star-formation is highlighted with a gray line for the median SFR of the sample (±0.4 dex scatter from Mineo et al. 2014), and may dominate the total X-ray observed emission for the X-undetected starbursts. Center panel: X-ray attenuation LX, obs/LX, int as a function of the infrared-based attenuation AV, tot (in a mixed model geometry and toward the center) for our sample of mid-IR detected AGNs. We assumed here a bolometric correction factor LX, intr, AGN = 0.04 × LBOL, AGN (Vasudevan & Fabian 2007). Stacks on the whole sample and on the X-undetected subset are displayed with hatching circles, while the violet shaded regions indicate the area of no obscuration, which incorporates a factor of two uncertainty in the conversion between intrinsic X-ray and bolometric AGN luminosity. Right panel: same diagram as before, but assuming an LBOL, AGN- dependent bolometric correction (Lusso et al. 2012), as explained in the text.

In the text

	Fig. 12. LBOL, AGN/M* vs. distance from the mass-size relation of LTGs (distLTG) for our SBs sample. The Eddington limit is shown with a blue horizontal line, while the shaded area takes into account the spread of M* among our sample and the uncertainty of the relation between stellar mass and BH mass by Reines & Volonteri (2015). The Eddington ratio is shown on the right y-axis.
In the text

Fig. B.1. Diagram of [O III]/Hβ compared to [S II]6717+6731/Hα (top) and [O III]5007/[O II]3727+3729 vs. [N II]6583/Hα (bottom) for 8 galaxies of our sample with optical spectra available and contemporary detection of [O III] and one between [S II]+Hβ or [O II]. The two diagrams are also called the S2 and the O2 BPT diagrams, respectively. As in the classical BPT diagram in Fig. 5, they show that our starbursts have different line excitation properties, and those with higher line velocity widths are generally shifted toward the AGN and LINER regions, according to empirical separation lines derived in the local (z ≲ 0.3) Universe.

In the text

Fig. B.2. SED fitting of the 25 starburst galaxies studied in this paper, presented in numerical order. The SED fitting procedure is described in Sect. 2.1. Blue and red curves show the stellar component (Bruzual & Charlot 2003) and the AGN torus emission (Mullaney et al. 2011), while the dust continuum emission is shown in green (Magdis et al. 2012). Downward arrows show the 2σ upper limit photometry at given wavelength. In each panel, zspec is the spectroscopic redshift from our Magellan spectra, zphot is the photometric redshift from Laigle et al. (2016), z is the output redshift of the best-fit SED (which in our cases was fixed to the spectroscopic value), and SFR is the best-fit infrared star-formation rate. The vertical dotted lines indicates the optimized range for fitting the stellar SED only (in order to avoid contamination from AGN torus and dust emission). The rightmost line also indicates the starting wavelength for fitting the AGN + dust SED in our paper. We remind the reader that the radio photometry (at 1.4 and 3 GHz) was not used in the fit, as explained in Sect. 2.1.

In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

Fig. B.3. Comparison between the EWs of Hδ, Hβ, Hα and Paβ lines. The blue lines represent the best-fit linear relation, while the blue shaded areas indicate the ±1σ dispersion of the galaxies around the best-fit line. We display in top-left corner of each panel the Spearman correlation coefficient “R” and corresponding “p”-value, the residual χ2 of the correlation (“${\chi }_{\text{red}}^2$”) and the 1σ scatter around the best-fit line, whose equation (comprising 1σ errors on the best-fit coefficients) is shown in the bottom-right corner. In the four diagrams, all the points with available measurements (excluding upper limits) were used in the fit.

In the text

Fig. C.1. For each galaxy, in two contiguous panels (top and bottom): Hα+[N II]λ6583 Å and Paβ emission lines from the Magellan-FIRE spectrum (Paγ is shown for the galaxy ID 245158, since this line is used to infer the attenuation). Each panel is further divided into three sections: in the upper part is shown with a black line the 1D portion of the spectrum (fully calibrated) close to Hα or Paβ line, along with the best-fit Gaussian superimposed (red line). In the second part below we display the noise spectrum for the same spectral range of the object (green line), while the fit residuals (noise normalized, as described in Calabrò et al. 2018) are shown in the bottom section with a black line. The galaxies are presented in numerical order for double detections, while those with only Hα+[N II] detection are shown at the end.

In the text

	Fig. C.1. continued.
In the text

	Fig. C.1. continued.
In the text

	Fig. C.1. continued.
In the text