ALMA Reveals a Stable Rotating Gas Disk in a Paradoxical Low-mass, Ultradusty Galaxy at z = 4.274

With the latest optical catalogs from the HFF-DeepSpace group (Shipley et al. 2018), the upper limits from Herschel/SPIRE (Rawle et al. 2016), and the measured submillimeter/millimeter fluxes (Table 1), we correct the photometry for the known magnification (average of μみゅー = 7 over Az9 ¹⁰ ) and use MAGPHYS_highz (v2; Battisti et al. 2020) to model the full spectral energy distribution (SED). MAGPHYS (da Cunha et al. 2008, 2015) is designed to self-consistently determine galaxy properties based on an energy balance approach using rest-frame UV through radio photometry in a Bayesian formalism. In brief, MAGPHYS uses the stellar population models of Bruzual & Charlot (2003), assumes a Chabrier (2003) initial mass function, and uses the dust model of Charlot & Fall (2000). We refer readers to the references above for more details. Simulations have shown that MAGPHYS does a good job of recovering the physical parameters of isolated galaxies and major mergers except during the near-coalescence phase (Hayward & Smith 2015).

The best-fit SED model is shown in Figure 2, and the best-fit parameters and their uncertainties are given in Table 1. It is reasonable to question the energy balance assumption in the SED fitting, especially since we find spatial offsets between the optical and infrared light (top left panel of Figure 1; see discussion in Section 4.3). In order to test this, we refit the SED excluding the IR data points. All parameters derived from the SED fitting, including the stellar mass, are completely consistent with the fits that include the IR points (see distributions in the bottom panel of Figure 2). Interestingly, we find that when the IR bands are excluded, the best-fit SED still lines up perfectly with the ALMA points and predicts the same IR luminosity. This might be surprising, but it seems that the bands between the Lyman limit and Lyαあるふぁ line (F555W, F606W, and F625W) are setting strong constraints on the dust. This is shown as the orange curve, which is the best fit excluding the ALMA and F555W, F606W, and F625W data points. This appears to be a consequence of the dust attenuation curve parameterization used in MAGPHYS and supports the energy balance argument even though there are offsets in the emission regions.

Zoom In Zoom Out Reset image size

Outputs from MAGPHYS include the stellar mass, star formation rate (SFR_SED), and IR and UV luminosities (L_IR, L_FUV). The SFR_SED is the sum of the stellar mass formed in the last 100 Myr and can be compared to the sum of the unobscured and obscured SFRs estimated from L_FUV and L_IR, respectively.

The L_FUV is calculated by fitting the UV continuum (1250 Å ≤ λらむだ ≤ 2600 Å) as a power law with βべーた as the UV slope. The UV slope based directly on the photometry is βべーた_phot = −1.10 ± 0.26. We obtain L₁₆₀₀ = (3.09 ± 0.36) × 10¹⁰ L_⊙ and convert this to FUV based on L_FUV ∼ 0.97 × L₁₆₀₀. The L_IR(8–1000 μみゅーm) is calculated by converting L_dust from MAGPHYS, which is the integrated dust emission at all wavelengths. The SFR_UV and SFR_IR are then calculated using the relations found in Murphy et al. (2011). All values are given in Table 1.

The sum of the obscured and unobscured SFRs, SFR_IR + SFR_UV = 30.3 M_⊙ yr⁻¹, is only 15% larger than SFR_SED. This can be attributed to the different assumptions inherent to each method (e.g., Utomo et al. 2014) and suggests that the energy balance assumption in MAGPHYS is reasonable.

We use the public lensing models provided through the HST Frontier Fields program ¹¹ to perform source plane reconstructions. Specifically, we use the CATS noncored model (Limousin et al. 2016), which predicts a redshift of 4.1 ± 0.2 for Az9, consistent with our new spectroscopic redshift from [C ii]. In HST, Az9 is extended north–south in the image plane over about 3'', and the dust and [C ii] extend another 2'' to the south. The range in the magnification across the source (in HST and ALMA) is μみゅー = 6–8.5, and the average magnification is μみゅー = 7 ± 1. The region used for the reconstruction of the HST image is shown with red dashed lines in the top left panel of Figure 1; it includes all of the clumpy faint emission that may be associated with Az9 and excludes another known multiply imaged source, 12.2.

Source plane reconstructions at a redshift of z = 4.27 are performed with Lenstool (Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009) using the publicly available CATS parameter files, which explicitly contain the optimized set of lens mass profiles (also determined with Lenstool). For each image, the pixels in the image plane are oversampled by a factor of 16 to provide greater accuracy, and source plane pixels are chosen to oversample by a factor of 8 (reflecting the magnification provided by lensing). For spectral cubes, the reconstruction is performed channel by channel.

We generate reconstructed images of [C ii], the Band 7 continuum, the Band 6 continuum, and the HST H-band image (bottom panels of Figure 1). In addition, we reconstruct the beam for each ALMA image by placing the observed synthesized beam in the center of the lensed image. As the magnification varies only modestly across the extent of the source, the source plane beam likewise does not vary much in the region of interest.

We measure the continuum fluxes on the reconstructed source plane images and confirm that the integrated fluxes are consistent with the integrated flux in the image plane maps scaled by the average magnification.

The bottom right panel of Figure 1 shows that the velocity map of [C ii] has a clear rainbow pattern consistent with rotation. We quantify this by performing a kinematic analysis on the source plane reconstructed cube using the ^3DBAROLO software (BBarolo; Di Teodoro & Fraternali 2015), which fits tilted ring models to emission-line data cubes. BBarolo models the geometric properties of the galaxy, which include the geometric center (x₀, y₀), inclination (i), and position angle (PA). The kinematic properties are the maximum velocity (V_max) and dispersion (σしぐま₀) of each ring.

First, we run the SEARCH algorithm in BBarolo, which is based on the DUCHAMP algorithm (Whiting 2012), to create a noise map. This is equivalent to masking pixels less than three times the rms in the moment 0 map. We limit the maximum ring radius to be within the SEARCH noise map, approximately 035 along the major kinematic axis. Prior to fitting the data, we assume a thin disk by fixing the scale height of the rings to z^o = 001; the resulting kinematic parameters are unchanged if we assume a 2×, 2.5×, or 10× thicker disk. The initial guesses for V_max and σしぐま₀ are 100 and 50 km s⁻¹ based on an initial analysis of the moment 1 and 2 maps, but we note that the final derived quantities are not sensitive to these choices. We use a two-stage fit to the data by first fitting both the kinematic and geometric parameters using a ring width equal to the spatial resolution, r_ring = 01, along the major kinematic axis with a reconstructed point-spread function (PSF) FWHM = 024. This ensures a good fit to the geometric center, inclination, and PA given the spatial resolution of our reconstructed image (Su et al. 2022). The best-fit inclination angle and PA are 466 and 1892, respectively. We take the best-fit parameters from this fit as inputs into a subsequent fit, where we only allow the maximum velocity and dispersion to vary and decrease the ring width by a factor of 2 to better sample the rotation curve. Throughout these fits, we approximate the PSF as a 2D Gaussian with a major/minor FWHM and rotation angle derived from a fit to the reconstructed source plane PSF, which is then convolved with the rotation model prior to calculating the residuals. As noted in Section 3.2, the source plane PSF varies little over the extent of the line emission because of a modest change in magnification.

The best-fit model is shown in Figure 3, and the derived kinematic parameters are listed in Table 1. The half-light radius (r_1/2) is calculated from the azimuthally averaged radial profile of the moment 0 model with errors propagated from the data and the fit. Here V_max is the average rotational velocity of all rings with $r\gt 0\buildrel{\,\prime}\over{.} 05$ because the inclination-corrected rotation curve is consistent with being flat at these radii. We do not expect high gas dispersion away from the core of the galaxy if it is disk-dominated, and indeed, the gas dispersion in the most extended rings is zero within their respective uncertainties. Therefore, we take the average dispersion over all rings with nonzero dispersion (all at $r\lt {r}_{1/2}=0\buildrel{\,\prime}\over{.} 26$) to be σしぐま₀. We estimate 1σしぐま uncertainties on the derived quantities following the Monte Carlo scheme built into BBarolo, which randomly samples models around the minimized residual (see Di Teodoro & Fraternali 2015 for further details). For parameters like V_max and σしぐま₀, we propagate the uncertainty per ring into the final averaged quantity's error.

Zoom In Zoom Out Reset image size

The moment 2 residuals in Figure 3 show a peak in the northern half of the galaxy, possibly arising from outflowing gas not associated with disk rotation. We refit the kinematic parameters using only the southern/redshifted portion of the galaxy and find that the region of high dispersion persists. The dispersion residual is a factor of 2 lower than when fitting the entire moment map and within the measurement uncertainty. Given that the residuals improve when masking this region, the dispersion peak in the moment 2 map in Figure 3 is most likely associated with measurement uncertainty. We conclude that our data do not have a sufficient S/N to support the presence of outflowing gas.

The bright [C ii] detection provides a spectroscopic redshift of 4.2738 for Az9 (top right panel of Figure 1). Given the brightness of the line and the redshift priors, this line can only be identified as [C ii]. The [C ii] spectroscopic redshift is consistent with the previous redshift estimates from optical photometry (Pope et al. 2017) and lens modeling (Diego et al. 2015; Limousin et al. 2016). Previous optical spectroscopic efforts failed to identify a redshift for this UV-selected galaxy. Therefore, [C ii] may be one of the best ways to identify redshifts in lower-mass dusty galaxies.

As discussed in Pope et al. (2017), the AzTEC beam (85) covered both source 5.2 (also known as Az9) and another multiply imaged galaxy, source 12.2, which is at z_spec = 1.71. With the high spatial resolution ALMA detections, we can now definitively rule out any millimeter emission coming from 12.2 (top left panel of Figure 1). Furthermore, the line we attribute to [C ii] cannot be from 12.2, since there are no known lines at the corresponding rest frequency.

Treu et al. (2015) reported a probable HST grism redshift of 0.928 for 5.1, the second-brightest image of this system. This redshift is inconsistent with the line we detect in the ALMA spectrum for 5.2, and we conclude that the grism redshift for 5.1 is incorrect.

From the MAGPHYS SED model (Figure 2), the intrinsic best-fit stellar mass is log(${M}_{\mathrm{star}})={9.33}_{-0.01}^{+0.21}$, sitting below the estimated knee in the stellar mass function at z ∼ 4 (Muzzin et al. 2013) and probing an unexplored region of stellar mass at this epoch. The stellar mass determined in this work, with the spectroscopic redshift and updated photometry, is slightly lower, although consistent within the uncertainties, than the stellar mass calculated using a different SED fitting code in Pope et al. (2017). With an sSFR = 12 Gyr⁻¹, Az9 is within the scatter of the star formation main sequence for its redshift (left panel of Figure 4).

Zoom In Zoom Out Reset image size

The width of the [C ii] line (282 ± 18 km s⁻¹) is consistent with the average found for UV-selected galaxies at z = 4–6 from ALPINE (Béthermin et al. 2020) and much narrower than the millimeter-selected galaxies from SPT (Gullberg et al. 2015). We measure a total intrinsic [C ii] luminosity of 4.63 ± 0.18 × 10⁸ L_⊙ and L_{[C II]}/L_IR = 0.0027. Az9 is not deficient in [C ii] like similarly dusty high-redshift galaxies (e.g., Gullberg et al. 2015) and has L_{[C ii]}/L_IR consistent with measurements of star-forming regions in nearby star-forming galaxies and the high-redshift ALPINE sample (Smith et al. 2017; Schaerer et al. 2020).

For its redshift, Az9 is a typical galaxy in terms of L_{[C ii]}, SFR, and stellar mass. However, unlike most low-mass galaxies, Az9 is very dust-obscured. The right panel of Figure 4 shows the fraction of the SFR that is obscured by dust as a function of stellar mass. Az9 is above the best-fit relation for z = 0–2.5 (green), and its f_obscured is >4× higher than the relation fit to the full ALPINE sample at z = 4–6 (blue dashed curve fit to open blue triangles; data from Béthermin et al. 2020). Put another way, the stellar mass of Az9 would need to be an order of magnitude larger to sit on the f_obscured relation for the ALPINE sample. Additional systematic uncertainties from SED fitting and lens modeling are insufficient to underestimate the stellar mass so severely; Az9 would remain overly dusty for its stellar mass.

Az9 is clearly detected in the Band 6 and 7 continuum images (left panels of Figure 1). In both cases, the emission is resolved, and we measure the fluxes through an optimal aperture in the image plane. The observed fluxes are presented in Table 1.

In both the image and source plane, the dust continuum emission is offset from the HST H-band emission, suggesting that the dust-obscured activity and stellar emission are coming from different regions of the galaxy. The HST emission is probing a rest-frame wavelength of 0.3 μみゅーm, which is highly susceptible to dust attenuation and may give rise to spatial offsets. While there is some faint HST emission under the [C ii] and dust continuum centroids and a bit of clumpiness to the north, the bulk of the HST emission is centered on a region offset by 02 (1.4 kpc) from the center of the dust continuum contours. This offset is consistent with the average offset found between UV and IR emission in galaxies in ALPINE (Fujimoto et al. 2020). Interestingly, the unobscured stellar light is coincident with the blueshifted [C ii] emission (bottom right panel of Figure 1). While it is possible that an outflow is clearing the dust on the blueshifted side of the disk, we see no evidence of this in our kinematic analysis (Section 3.3).

From the kinematic modeling, we obtain a [C ii] half-light radius of 026, which corresponds to 1.8 kpc at z = 4.274. This radius is consistent with the range of [C ii] sizes measured for lower-mass galaxies in ALPINE (Fujimoto et al. 2020).

Figure 3 shows the results of our kinematic modeling. Az9 shows a smooth velocity gradient across the galaxy, defining the kinematic axis, and a centrally peaked velocity dispersion distribution. We calculate V/σしぐま = 5.3 ± 3.6, which is consistent with a stable, dynamically cold, rotating disk. We compare our measured V/σしぐま with values found in simulations and other observations but note the caveat that different definitions and ways of measuring V and σしぐま may lead to different values. Az9 has a V/σしぐま that is >2× higher than the V/σしぐま predicted for similarly low-mass galaxies at z = 4.3 from the TNG50 simulation, driven primarily by the higher velocity dispersion predicted in the simulation (Pillepich et al. 2019). With clear rotational disk kinematics, Az9 does not show any indications of being a merger.

Rizzo et al. (2022) tested the robustness of V/σしぐま at classifying disks and mergers using mock [C ii] observations and found that there are some cases where a merger can look like a disk in V/σしぐま; however, the mergers are always picked out with multiple-peaked emission profiles. For our kinematic analysis, we have ∼1.5× resolution elements over the major kinematic axis; in this regime, disk galaxies are always classified correctly, but ∼50% of mergers are misclassified as disks (Rizzo et al. 2022). While we cannot rule out a merger based on V/σしぐま alone, the [C ii] profile of Az9 is nicely single-peaked (Figure 1) in favor of the disk interpretation.

There are only handfuls of rotating disks observed at z > 4 and even fewer at low stellar masses. Observations show that the fraction of rotation-dominated (disk) star-forming galaxies with a stellar mass of <10¹⁰ M_⊙ at z = 3 is <40% (Förster Schreiber & Wuyts 2020), with the disk fraction increasing with stellar mass (Tiley et al. 2021). Of the 29 [C ii]-detected galaxies in ALPINE with detailed kinematic modeling with BBarolo, only 6 (21%) were classified as rotators with stellar masses of ∼10¹⁰ M_⊙ (Jones et al. 2021). Rizzo et al. (2020) presented a remarkable rotator at z = 4.2, surprisingly unaffected by a nearby companion (Peng et al. 2023) but with a stellar mass 6× higher than Az9 (see also Roman-Oliveira et al.2023). Isobe et al. (2022) found that local galaxies with low masses (<10⁹ M_⊙) are all observed to have V/σしぐま < 1. Therefore, Az9 appears to be an outlier compared to existing observations with higher V/σしぐま for its stellar mass and redshift.

In addition to showing clear rotation, the low velocity dispersion of Az9 (26 km s⁻¹) suggests that it is stable. Dusty star-forming galaxies at z ∼ 2 from the KAOSS survey have average rotational velocities and velocity dispersions from ionized gas of 190 and 90 km s⁻¹, respectively (Birkin et al. 2023). Even though these dusty galaxies are technically rotation-dominated, their high dispersion suggests turbulent rotating disks. Even accounting for the fact that the ionized gas dispersion is observed to be ∼10–15 km s⁻¹ higher than the molecular/atomic gas dispersion in galaxies out to z ∼ 2.6 (Übler et al. 2019), Az9 still has a much lower velocity dispersion. This may be expected for its lower mass but perhaps unexpected given its high dust content. It is unclear what role dust plays in disk turbulence, and measuring the kinematics of gas in multiple phases in higher-redshift and lower-mass galaxies will help address this question.

From the rotational velocity and the radius, we calculate a dynamical mass at 2 × r_1/2 of M_dyn = 1.6 × 10¹⁰ M_⊙, which gives M_*/M_dyn ∼ 0.1. Measurements of the molecular gas mass are needed to complete the census of baryonic mass in Az9 and constrain the dark matter fraction.