© The Authors 2025

1. Introduction

Galaxies at redshifts of z > 6 represent the early stages of galaxy formation and evolution, often referred to as cosmic dawn. Discovering these galaxies and determining their physical properties is key to understanding the further phases of galaxy evolution (e.g., Bouwens et al. 2022a; Finkelstein et al. 2023; Rieke et al. 2023). Furthermore, these galaxies are essential for understanding the galaxy-scale processes that gave rise to the epoch of reionization (e.g., Hutter et al. 2021). The gradual change of a steeper slope of the faint end of the ultraviolet (UV) luminosity function with increasing redshift (e.g., McLure et al. 2013; Finkelstein et al. 2015; Oesch et al. 2018; Bouwens et al. 2022b) implies an increased relative importance of the sub-L* population’s contribution to the overall production of the UV light responsible for reionization. However, mapping galaxies at redshifts z > 6 is challenged by the enormous distances resulting in faint appearance, the redshifting of the rest-frame UV/optical/near-infrared (NIR) light into wavelengths with low atmospheric transmission, and the varying morphologies (i.e., clumpiness; e.g., Gullberg et al. 2018) of galaxies in the early Universe might be more difficult to resolve with sensitive interferometers.

An efficient way to study the sub-L* z > 6 galaxy population is through gravitational lens studies. The gravitational lensing magnification aids the time expensive observations of the low-luminosity galaxy population enabling detailed studies of individual sources. Using submillimeter wavelength observations to resolve emission lines is a necessary means for determining physical properties and characterizing the dynamical state of such galaxies. The far-infrared (FIR) fine-structure lines of ionized carbon and oxygen has proven essential for this (e.g., Hashimoto et al. 2019; Harikane et al. 2020). For local galaxies, even dwarf galaxies, relations between the star formation rate (SFR) and the line luminosity of the FIR fine-structures lines have been established (e.g., De Looze et al. 2014; Cormier et al. 2015). However, the growing number of [C II] observations of z > 6 sub-L* galaxies have yielded a more complicated result for the early epoch (e.g., Capak et al. 2015; Knudsen et al. 2016; Bradač et al. 2017; Carniani et al. 2017; Hashimoto et al. 2018; Tamura et al. 2019; Laporte et al. 2019; Bakx et al. 2020; Jolly et al. 2021; Fujimoto et al. 2021; Molyneux et al. 2022); for example, some of those early [C II] observations yielded several non-detections or line luminosities lower than initially predicted. With the combination of high sensitivity and sub-arcsecond resolution of modern submillimeter interferometric observatories, it has become possible to resolve line emission even for z > 6 galaxies. This makes it possible to conduct dynamical studies (e.g., Carniani et al. 2013; Litke et al. 2019; Spilker et al. 2022; Rowland et al. 2024); however, without gravitational lensing, reaching sub-kpc scales remains challenging. Numerical zoom-in simulations of z > 6 galaxies show that the kinematics can be very complex. For example, Kohandel et al. (2019) demonstrate that the same galaxy over the course of hundred million years can change appearance from a merging system to a rotating disk. The same simulations, as well as those of Katz et al. (2022), also reveal how the predicted observations are dependent on the orientation, showing the difficulty of distinguishing different dynamical states with only low angular resolution.

The galaxy A1689-zD1 was reported as a z > 7 candidate with a photometric redshift of z = 7.6 (Bradley et al. 2008) and this was confirmed with X-shooter spectroscopy that showed a break in the continuum emission at a redshift of z = 7.5 ± 0.2 (Watson et al. 2015). A1689-zD1 was also the first sub-L* with a dust-detection with the Atacama Large Millimeter/submillimeter Array (ALMA), revealing not only a dust mass similar to that of the Milky Way interstellar medium (ISM), but also a complex distribution that could be interpreted as a potentially interacting system (Watson et al. 2015; Knudsen et al. 2017).

Some of the brightest FIR fine-structure lines are [C II] 158 μm and [O III] 88 μm in star-forming galaxies. As carbon has an ionization potential of 11.2 eV, which is below that of hydrogen, the [C II] line is used for both tracing neutral and ionized gas. Ionized oxygen has an ionization potential of 35.1 eV and, thus, it is only observed in the ionized medium (e.g., H II regions). Thus, while each of these two bright emission lines can be used to resolve the properties of the ISM, when used in combination, they also provide insights to the ionization state of the ISM in early galaxies (e.g., Harikane et al. 2020; Killi et al. 2023).

Here, we present deep imaging of the two FIR fine-structure lines [C II] 158 μm and [O III] 88 μm for A1689-zD1. Some of these data were recently analyzed and the results have begun shedding light on different aspects about FIR properties (Bakx et al. 2021; Akins et al. 2022; Wong et al. 2022; Killi et al. 2023). While the original goal of the [C II] observations was to determine the redshift, the [C II] and [O III] data enable imaging of sub-kpc scales thanks to the gravitational lensing amplification (thereby enabling multi-component decomposition). Furthermore, we present detection limits from molecular gas observations. The paper is structured as follows. Section 2 presents the ALMA observations and data, along with Green Bank Telescope (GBT) observations of CO(3–2). Section 3 gives the results including an analysis of the complex gas kinematics. A discussion is given in Sect. 4 and our conclusions are summarized in Sect. 5. We assume a Λ cold dark matter cosmology, with H₀ = 67.7 km s⁻¹ Mpc⁻¹, Ω_M = 0.307, and Ω_Λ = 0.693 (Planck Collaboration XIII 2016).

2. Observations

2.1. ALMA data

Observations of A1689-zD1 were obtained with ALMA to search for the [C II] 158 μm and subsequently to follow-up the [O III] 88 μm line (ALMA projects 2015.1.01406.S and 2017.1.00775.S). The [C II] observations were tuned to match frequencies of candidate emission lines found in previous data (Knudsen et al. 2017). The receivers were tuned to 217.19–220.94 and 231.98–235.74 GHz, where the line was predicted to be in the former frequency range. The calibrators used were J1337−1257 for flux calibration, J1256−0547 for bandpass, and J1319−0049 for gain calibration. Data were calibrated using the ALMA pipeline; we verified the observatory delivered calibration with only the need for limited additional flagging, however, three out of ten executions needed the absolute flux calibration redone as the pipeline had not included the correct values for the flux calibrator. After this correction, the uncertainty of the absolute flux calibration is estimated to be about 5%¹. Combining the data from all ten executions, we reached an overall depth of 0.18–0.2 mJy/beam in a 5 km s⁻¹ channel² when imaging using a Briggs weighting with a robust parameter of 0.5; the angular resolution is 0.29″ × 0.27″ and the position angle (PA) is −80.8°. When combining the lower frequency band with the continuum data from Knudsen et al. (2017), we obtained a root mean square (rms) flux uncertainty of about 6.3 μJy/beam.

The emission line in the upper frequency band was detected at ν = 233.3 GHz and identified as the [C II] 158 μm emission line corresponding to a redshift of z_[CII] = 7.133. Using this redshift ALMA band 8 observations were obtained for the [O III] 88 μm line. The receivers were tuned to 404.40–408.21 and 416.33–420.09 GHz. Calibration was performed using J1229+0203 for flux and bandpass calibration, and J1256−0547 for gain calibration. The observatory-delivered calibrated data were checked and considered sufficient, and no further treatment was done. The data were imaged, yielding a continuum map with rms of 29 μJy/beam and an line cube with an rms of 0.39–0.44 mJy/beam in ∼11 km s⁻¹ (15.625 MHz) channels. The resulting images have an angular resolution of 0.33″ × 0.28″, PA = $-87.3°$, which is similar to the resolution of [C II] and band 6 maps. A conservative estimate on the absolute flux calibration is 10%. A continuum subtraction was performed for both the [C II] and [O III] data cubes using the CASA task UVCONTSUB where the continuum was estimated using the line-free channels.

ALMA band 3 observations of the CO(6–5) and CO(7–6) lines were obtained using a receiver setup of 84.07–87.83 GHz and 96.16–99.94 GHz (ALMA project 2017.1.00775.S). The calibrator J1256−0547 was used for flux and bandpass calibration, while J1319−0049 was used for gain calibration. The observatory-delivered, calibrated data were of sufficient quality. The data were imaged using natural weighting and also using uv-tapering with a circular Gaussian profile of 0.5 arcsec resulting in an angular resolution of 0.92″ × 0.80″ (PA = 64°). The rms in the continuum map is 11 μJy/beam, and in the spectral cube near the frequency of redshifted CO(6–5) [CO(7–6)] the rms is 0.23 [0.20] mJy/beam in a 15.625 MHz channel (corresponding to a 55 [47] km s⁻¹ channel width).

2.2. GBT data

Using the GBT³ observations, we carried out a search for the redshifted CO(3–2) emission with the Q-band (40–50 GHz) receiver. Six observing sessions (AGBT14A_565) were undertaken in September-October 2014 and an additional eight sessions (AGBT18A_230) were carried out in February–September 2018. In total, 64 hours of telescope time were used for searching for CO(3–2) emission. The data were collected with the VEGAS spectrometer using the wide bandwidth modes 1 and 3. The observations were done using the sub-reflector nodding (“SubBeamNod”) observing mode with a six-second switching period between the two Q-band beams that are separated by 58 arcsec on the sky. At the central observed frequency of 42.5 GHz, the full width at half maximum (FWHM) resolution of the GBT beam is about 17 arcsec. The bright nearby quasar 3C279 was used for pointing and focus, while 3C286 was used for absolute flux calibration (Perley & Butler 2017). After correcting for atmosphere, pointing drifts, and calibration, the uncertainty on the final flux density scale is approximately 15%.

2.3. HST data

We reprocessed all available Hubble Space Telescope (HST) and Spitzer data in archive in the A1689 field as a part of the Complete Hubble Archive for Galaxy Evolution (CHArGE; Brammer, in prep.), which aims to perform uniform processing and analysis of all archival HST and Spitzer data taken away from the Galactic midplane (Fujimoto et al. 2022). We aligned all of the HST exposures to sources in the Gaia DR2 catalog (Gaia Collaboration 2018) and created final mosaics in a common pixel frame with 50 mas and 100 mas pixels for the Advanced Camera for Surveys (ACS) Wide Field Channel (WFC) and Wide Field Camera 3 (WFC3) IR filters, respectively. We aligned the individual Spitzer exposures to the same astrometric frame as the HST frame and generated final drizzled IRAC mosaics with $0\stackrel{″}{.}5$ pixels. Further details of the HST (Spitzer) image processing with the grizli (golfir) software is presented in Kokorev et al. (2022).

3. Results

3.1. [C II] and [O III]

Emission lines are detected in both ALMA data sets consistent with [C II] 158 μm and [O III] 88 μm. Both lines show a complex line profile as well as complex surface brightness distribution. Similarly, the continuum emission shows a complex structure consistent with the integrated emission line distribution. In Fig. 1, we show the moment-0 and moment-1 maps, while in Figs. 2 and 3, we show the line emission per 40 km s⁻¹ channel for both lines. The integrated line intensity along with continuum flux densities are given in Table 1. The band 6 and band 8 continuum results are consistent with the photometric results from Knudsen et al. (2017) and Inoue et al. (2020).

Fig. 1. Velocity-integrated maps of [C II] (left) and [O III] (right) of A1689-zD1 with the moment-0 (top) and moment-1 (bottom), derived from the velocity ranges v[CII] = −215 − 270 km s−1 and v[OIII] = −250 − 270 km s−1, respectively. The contours represent the 3, 5, 8, 12, and 17σ levels, where σ[CII] = 11.0 mJy km s−1 beam−1 and σ[OIII] = 63.2 mJy km s−1 beam−1. The moment-1 map is masked at 3σ of the moment-0 map.

Fig. 2. Mosaic of the [C II] binned by velocity in 40 km s−1 intervals. The value in the upper right corner of each sub-image indicates the middle of the velocity relative to z = 7.133. The gray contours show the band 6 continuum, where the levels correspond to 4, 6, 8σ. The five components of the analysis (Sect. 3.5) are indicated by name. The color bar below the last panel shows the level corresponding to the velocity integrated flux, SdV, for each map.

Fig. 3. Mosaic of the [O III] observations. The upper right corner shows the central velocity of each 40 km s−1 channel relative to z = 7.133. The gray contours show the band 8 continuum, where the levels correspond to 4, 6, 8, 10, 12, and 14σ. The five components of the analysis (Sect. 3.5) are indicated by name. The color bar below the last panel shows the level corresponding to the velocity integrated flux, SdV, for each map.

3.2. Redshift

The [C II] 158 μm and [O III] 88 μm lines suggest central redshifts of 7.1329 and 7.1334, respectively, based on a single Gaussian profile fit. Taking the mean of these as the systemic redshift and the difference as the uncertainty, we find A1689-zD1 is at $z=7.1332\pm 0.0005$. This is significantly (∼13 000 km s⁻¹) offset from the redshift derived from the NIR spectrum of the Lyman-break in Watson et al. (2015), which was 7.5 ± 0.2. This is similar to the 11 000 km s⁻¹ difference seen in redshift determination using optical/NIR spectroscopy and ALMA spectroscopy for MACSJ1149-JD at z = 9.11 in [O III] emission (Hashimoto et al. 2018) and the photometric redshift of 9.51 (Hoag et al. 2018), as well as for GNz11 of z_grism = 11.1 versus z_spec = 10.6 confirmed with the JWST (Oesch et al. 2016; Bunker et al. 2023). As bright Lyman-alpha emission or other UV emission lines have not been detected in the spectra, the optical redshift has been derived based on the continuum break. However, given the potential large column of HI along the line of sight as well as the relative faintness of these distant galaxies provides a challenge for obtaining redshifts. A detailed analysis of this will be provided in Watson et al. (in prep.).

3.3. Molecular gas

Figure 4 shows the result of the GBT observations. The CO(3–2) line was not detected, and the 3σ limit is given in Table 1 assuming a linewidth of 300 km s⁻¹. To mitigate residual baseline structures high-pass filtering methods were applied to the individual SubBeamNod data scans. These filtering methods would remove any underlying broad features (> 1000 km s⁻¹), while preserving emission for the expected line width of about 300 km s⁻¹ and, thus, resulting in a flat spectrum (Fig. 4). Although the line was not detected, if we were to smooth the data further, there is a weak, nearly 2σ bump (with a width of about 200 ± 100 km s⁻¹ at z = 7.13 when fit by a Gaussian) at the expected frequency. Hundreds of hours of additional GBT time would be required to test whether this “bump” is real.

Fig. 4. Non-detection of CO(3–2) emission measured at the GBT smoothed to a channel resolution of about 30 km s−1. The expected line location based on the [C II] and [O III] redshift of z = 7.133 is shown by the vertical dotted line.

The CO(6–5), CO(7–6), and the [C I] (1–0) 492 GHz lines were not detected in the ALMA band 3 data. Assuming a linewidth of 300 km s⁻¹ and an area similar to the [O III] detection, we estimated 3σ upper limits for all three transitions. These are given in Table 1. We note that detection of moderately high-J CO transitions is a challenge of z > 6 galaxies (e.g., Hashimoto et al. 2023).

3.4. Gravitational lensing

We used a parametric model of the Abell 1689 mass distribution to estimate the magnification and distortions at the location of A1689-zD1, as well as across the source. Our reference mass model is constructed with the Lenstool software (Jullo et al. 2007) and described in Limousin et al. (2007).

At the measured spectroscopic redshift of z = 7.133 this model predicts a single image at this location with a total magnification of 9.6, involving an angular magnification of ∼6×, along a specific direction at PA = 85 degrees (roughly east-to-west), and 1.6× perpendicular to it. The amplification factor and shear directions do not vary by more than 10% and 5°, respectively, across the 1.5″ extent of the source (Table 2). This is due to the fact that there is no significant cluster substructure or cluster member identified in the vicinity of A1689-zD1, which could locally affect the mass distribution and magnification.

To estimate the errors on the magnification factor, we used the best-fit mass model and the range of Markov chain Monte Carlo samples from the optimisation performed in Lenstool by Limousin et al. (2007). We obtain a statistical relative error of 19% (Table 2). In addition, we make use of additional existing mass models of the same cluster (e.g., Li et al. 2012) where a different value is used for the relative contribution of the cluster-scale and galaxy-scale distribution of the mass model, to estimate the systematic uncertainty on the average magnification of the source. We estimate a systematic relative error of 10%.

The angular magnification implies an angular resolution in the [C II] map of ∼0.037″ × 0.15″, which correspond to 0.2 × 0.78 kpc. Correcting the sensitivity for the magnification, this implies a spectroscopic rms of ∼20 μJy/beam per 5 km s⁻¹ channel and a continuum sensitivity of 0.67 μJy/beam for the band 6 data, making this some of the deepest data obtained for a single high-z galaxy with ALMA.

For the [O III] data, the magnification-corrected resolution is ∼0.051″ × 0.18″ corresponding to 0.27 × 0.95 kpc. The spectroscopic sensitivity corresponds to 42–47 μJy/beam in a ∼11 km s⁻¹ channel and the continuum sensitivity is 3.1 μJy/beam.

3.5. Decomposition into galaxy components

When stepping through the continuum-subtracted ALMA line datacubes, both [C II] (Fig. 2) and [O III] (Fig. 3) show emission peaking at different spatial locations in the source as a function of velocity. This complex structure is dominated by three main spatial components with a maximum flux at ∼ − 80 (“blue”), ∼0 (“green”), and ∼ + 160 km s⁻¹ (“red”) and present in both emission lines. This multi-component structure was already evident from the complexity seen in previous continuum results (Knudsen et al. 2017). However, emission lines at a high S/N now provide us with a better possibility for disentangling this scenario.

First, as made clear from the relative spatial locations of the main components, as well as the moment-1 map derived from the cube (Fig. 1), they do not follow a regular structure (as seen, e.g., in regular rotating disk galaxies). Therefore, we identified and model each component from the datacubes using an iterative process: (1) we produce a narrowband image from the line data comprising the main channels of emission in each component; (2) we fit this image with a 2D elliptical Gaussian model for the line flux distribution (including the central location, total flux, full width at half maximum (FWHM) along major and minor axes, and PA of the major axis); (3) we extracted the corresponding spectrum of this component and fit line parameters with a 1D Gaussian (central velocity, v₀, velocity dispersion, σ₀); and (4) we subtracted a 3D model of line emission for this component combining the spatial and spectral model.

Due to the proximity of some components (both spatially and in velocity space), this fit is performed three times iteratively, subtracting all other components from the datacube before a new fit. In addition to the three main components mentioned above, two fainter components peaking at ∼ + 80 (“yellow”) and ∼ − 220 (“violet”) km s⁻¹ were added, and more complex spatial distributions of the flux were used when necessary from the residuals. A few extended components also revealed a spatial variation of the line peak, which we fit using an empirical velocity gradient model (arctan model with a truncation velocity, v_t, and a truncation radius, r_t). We summarize the decomposition into five components below and we provide the full model parameters in Tables 3 and 4. We note that we originally identified the components from the [C II] data, but we recovered the same structure in the [O III] data, where we followed similar narrowband images and spatial distributions for the fit of each component.

Violet: The faintest of the five components. It was fit with a single non-rotating component and for [C II], the fit was made challenging by empty channels between the spectral windows.

Blue: Modeled with two elliptical Gaussian components. One is extended with rotation (denoted “blue1”) and one compact with only weak rotation (“blue2”). Their central positions are close but not identical.

Green: Broad component modeled with a single extended source, which has weak rotation.

Yellow: Single compact component with no measurable rotation.

Red: Combination of two individual elliptical Gaussian profiles with weak rotation, extended in the WSW-ENE direction.

Combining these five components, we recovered 95 and 97% of the total emission measured in the [C II] and [O III] line data, respectively, (Fig. 5) in the inner region (not the extended region studied in Akins et al. 2022). Some extra emission, somewhat diffuse, appears to be present but cannot easily be modeled as individual components. It is worth noting that the choice of decomposition (number of components and their parametrization), while somewhat arbitrary, is driven by the complexity and the high signal-to-noise ratio (S/N) of the data. Correcting the data for gravitational magnification, the extracted components are found to have sizes in the range 0.1–1 kpc.

Fig. 5. [C II] (top) and [O III] (bottom) spectra (continuum-subtracted) overlaid with the resulting models for each component (color coded according to their component names). The black line shows the sum of the models. The residuals after subtracting the model results are shown as a thin gray histogram.

For comparison of the derived model to the data, we show the moment-0 and moment-1 maps of both the observations and the model in Fig. 6 for both [C II] and [O III]. The residual maps indicate that the models effectively reproduce the observations. We note that for [C II] there could be a small offset and that this could be due to the [C II] emission having an additional extended component (e.g., Akins et al. 2022), which adds to the complexity of the source.

Fig. 6. Overall moment-0 maps (first and third rows) and moment-1 maps (second and fourth rows) for A1689-zD1 for [C II] and [O III] for the data (left), the model (middle), and the residuals (right). Contours show the model with the gray showing the moment-0 (3, 4, 5, 7, 10, 13, 16, and 19σ, where the rms is measured in the data regions excluding the source), and the white contours showing the moment-1 in steps of 20 km s−1. The applied mask is determined from the model. We note that this does not include the extended “halo” component presented in Akins et al. (2022).

In Fig. 7 we show the contours of the individual components from [C II] and [O III] overlaid on the HST images. The relative astrometry of the HST and ALMA images should be correct to $0\stackrel{″}{.}03$ allowing us to make a robust comparison of the galaxy’s UV and the ALMA results. As the HST data show no significant systematic color change across the galaxy, we used a stacked image combining the data from observations using the F125W, F140W, and F160W HST filters to increase the S/N of the rest-frame UV data. For further illustration of the data, we show the decomposed components in the 3D in Fig. 8, using an approach to showing 3D data similar to that of Spilker et al. (2022) and showing the emission of each component above S/N of 3.

Fig. 7. Contours of the individual components of the [C II] 158 μm (top) and [O III] 88 μm (bottom). Overlaid is the combined HST image obtained with F125W, F140W, and F160W filters. The contours are color-coded according to their component name; contours are in steps of 2σ starting at 4σ (see Fig. 14 for the rms).

Fig. 8. Decomposed components of the observations of A1689-zD1 shown in 3D with color-coding according to component name. The integrated moment-0 maps is shown at the end of the plot.

The clearest difference between the UV and FIR data is the number of components, with two components seen in the UV and the five components of the [C II] and [O III] data. Both the violet and yellow components have no obvious UV-counterpart. Given the apparent offsets between the green, blue, and red components, we cannot unambiguously associate the components with the UV-components. One possible interpretation is that green is associated with the SW UV-component and blue is associated with the NE UV-component. In contrast, the violet and yellow components do not show HST counterparts, while dust continuum emission is detected from these components. This suggests that completely dust-obscured regions exist in this system, consistent with the highly (≃80 − 90%) obscured SFR fraction of A1689-zD1 (Bakx et al. 2021; Akins et al. 2022). This is also reminiscent of the recent identification of completely dust-obscured HST-dark galaxies at z > 7 Fudamoto et al. (2021). We further discuss the possible origins of these components in Sect. 4 below.

3.6. Disentangling the dust continuum

We also estimated the flux densities of the continuum emission in band 6 and band 8 for each component. We first performed a forced-aperture photometry. Given that [O III] 88 μm shows a better correlation than [C II] 158 μm with the SFR (e.g., De Looze et al. 2014; Olsen et al. 2017), we optimized spatial positions and shapes of the apertures following the five components of the [O III] model. In the left panel of Figure 9, we show the five aperture positions and shapes and summarize the extracted flux densities in Table 5. In this approach, the apertures do not overlap with each other; however, we cannot rule out the possibility that some measurements within the aperture are affected by contributions from other components.

Fig. 9. Decomposition of the band 6 (about 158 μm, top) and band 8 (about 88 μm, bottom) dust continuum for the five components. The 4″ × 4″ cutouts of the observed, IMFIT model, and residual maps are presented from left to right. The solid and dashed contours show the intensity at the 2, 3, 4, 5, 6, 8, and 10σ levels and −3 and −2σ levels, respectively, which are measured based on the residual map. The spatial position of the five [O III] components are marked with the color crosses. The cyan circles indicate the elliptical apertures that are used for the aperture photometry.

Next, we modeled and decomposed each contribution by using CASA IMFIT. From the observed morphology, we assumed four components. The violet (red) component was not included in the fitting for the band 6 (band 8) maps, but we extracted its photometry with a radius of $0\stackrel{″}{.}2$ aperture on the residual map. We fixed the spatial positions of the faint components of red and violet when they are included in the model. We kept the flux density as a free parameter, while we fixed some of the other parameters (size, axis ratio, and PA) in a routine process to obtain stable results. In Table 5, we summarize the flux density estimates in the both bands for each component. By assuming a single modified blackbody and a spectral index of β = 1.7 in a spatially averaged manner (Akins et al. 2022), we also derived the dust temperature (T_dust) and IR luminosity (L_IR). Due to the low S/N for the violet and red components, we derived their L_IR values by fixing their T_dust = 41 K taken from the spatially integrated best-fit value presented in Akins et al. (2022). These estimates are also listed in Table 5. The sum of the five components gives 0.45 ± 0.20 mJy and 1.71 ± 0.58 mJy in bands 6 and 8, respectively. The band 8 continuum is comparable to the total flux estimate in Table 4, while the band 6 continuum model underestimates the total flux by ∼50%. This indicates that the diffuse continuum emission, significant in band 6 (Akins et al. 2022), is not fully recovered by this decomposition analysis. The continuum model and residuals are shown in Fig. 9. We note that we here use a higher resolution weighting scheme compared to Akins et al. (2022) and thus the difference is less evident in the residuals.

4. Discussion

4.1. Resolved and overall physical properties

Based on the component analysis from the previous subsections, we estimate the line luminosity and the continuum properties of each component. The line luminosity is estimated both using the observations and the model results. In Table 6 we provide the line luminosity derived for each component. The uncertainty is estimated from the observed data, but does not include the uncertainty on the absolute flux calibration.

Based on our resolved line and continuum estimates, we can evaluate physical properties of the full source and the individual components. While blue and red are further decomposed into two subcomponents in Sect. 3.5, here we focus on the physical properties of the brighter component that dominates their total emission.

4.1.1. Total SFR

We calculated the total SFR by summing the obscured (SFR_IR) and unobscured SFR (SFR_UV) based on the ALMA and HST observations. We measured the UV luminosity L_UV for blue and green components with a $0\stackrel{″}{.}2$-radius circular aperture placed at the peak positions of the two blobs observed in the HST maps. Because of the lack of the HST counterparts, we obtained a 3σ upper limit with the same size aperture for the other components. We used the same conversions from IR and UV luminosity to SFR as Akins et al. (2022) based on the calibration presented in Murphy et al. (2011) and obtain a total SFR of $16.5_{-6.8}^{+22}{M}_{\odot }$ yr⁻¹. We did not account for the calibration uncertainty in the error estimate.

We confirm consistency with previous measurements (Bakx et al. 2021; Akins et al. 2022) given the recovery fraction of the total dust continuum in our measurement (Sect. 4.1). We found that the obscured fraction f_obs is ∼80% in the blue and green components, respectively, and more than ∼90% in the other components, indicating the high obscuration in every part of the system. We summarize SFR_IR, SFR_UV, and f_obs for the individual components in Table 7.

4.1.2. SFR–L_line relation

In Fig. 10 we show SFR–L_[CII] and SFR–L_[OIII] relations for the individual components and the entire system of A1689-zD1. We also show for comparison results obtained in local galaxies and z > 6 galaxies without active galactic nuclei (AGNs) from the literature. For A1689-zD1, we find that the relation in both lines is consistent with previous results for z > 6 galaxies, where the [C II] line is slightly below the relation from local galaxies, but within its 1σ range; whereas the [O III] line agrees well with it. For individual components, we find that the relation has a little scatter, but still falls within the 1σ range of the local galaxies. These results are consistent also with ALMA results for z > 6 lensed galaxies (e.g., Knudsen et al. 2016; Fujimoto et al. 2021; Molyneux et al. 2022).

Fig. 10. Line luminosity vs. SFR for L[CII] (left) and L[OIII] (right). Black outline diamonds with color centers indicate the color-coded individual components from Sect. 3.5. The uncertainty on the SFR-axis for each component corresponds to the range given by a Tdust = 25 − 55 K. The orange and gray lines and associated regions show the relations for star-forming galaxies and low-metallicity dwarf galaxies from De Looze et al. (2014). Data from the literature is included for comparison for high-z galaxies without AGNs (Capak et al. 2015; Inoue et al. 2016; Knudsen et al. 2016; Pentericci et al. 2016; Bradač et al. 2017; Laporte et al. 2017, 2019, 2021; Carniani et al. 2017; Smit et al. 2018; Walter et al. 2018; Hashimoto et al. 2018, 2019; Tamura et al. 2019; Harikane et al. 2020; Binggeli et al. 2021; Fujimoto et al. 2021; Molyneux et al. 2022). We note that the data from the REBELS survey include the sources with L[CII] and UV + IR SED-estimated SFR (Sommovigo et al. 2022; Ferrara et al. 2022). On the right panel, for [O III] we also include the results from Marrone et al. (2018). In the left panel, the blue and magenta lines represent model predictions from Vallini et al. (2015) and Olsen et al. (2017).

Given that the emission is resolved, it is also natural to compare this to local galaxies. For example, the SHINING survey of nearby star-forming galaxies extensively studied FIR fine-structure lines and compared this to continuum properties (e.g., Herrera-Camus et al. 2018). Comparing line luminosity to the FIR luminosity and the FIR surface brightness, Σ_FIR, of A1689-zD1 to the results from the SHINING, we note that [C II] results are within the scatter and consistent with that of SHINING results, while the [O III] line appear about a factor two brighter than the brightest [O III] detections in the SHINING survey.

Fig. 14. Moment-0 and moment-1 maps for each component: [C II] results (left) and [O III] results (right). The models of the other components have been subtracted from the data cube. The left columns show the moment-0 maps integrated over the line of the component; the middle columns the moment-1 given relative to the A1689-zD1 overall velocity field with zero velocity corresponding to z = 7.133; and the right columns the moment-1 map for each individual velocity (same velocity used for both [C II] and [O III]). The overlaid contours represent the fit models, in the moment-0 maps these correspond to 3, 4, 5, 6, 7, 8, 9, and 10σ; in the moment-1 maps, these correspond to velocity in km s−1 in steps of 10 km s−1. The moment-0 maps are given in S/N, where the rms in each map is measured excluding the source region (for [C II]: rms = 5.4, 11.6, 11.5, 6.4, and 8.3 mJy km s−1 beam−1, respectively, and for [O III]: rms = 27.7, 50.8, 54.1, 29.4, and 38.1 mJy km s−1 beam−1, respectively (from top to bottom).

In a compilation of z = 5 − 7 galaxies with [C II] detections, Carniani et al. 2018 find that “normal” star-forming galaxies deviate from the local empirical relation of Σ_SFR − Σ_[CII] with a typically lower [C II] surface brightness. Comparing the SFR-L_[CII] relation from the ALPINE survey (Schaerer et al. 2020), assuming average extended emission of 8 kpc² (as done in Posses et al. 2025), it is positioned between the relation from Carniani et al. (2018) and the local relation from Herrera-Camus et al. (2015). Differences between the different populations (e.g., local vs. high-redshift galaxies) could be caused by heating mechanisms, effects of metallicity, or the ionization parameter (e.g., Carniani et al. 2018). In Fig. 11, we compare the resolved [C II] emission of A1689-zD1 with those relations and with recent results of the two resolved high-z galaxies REBELS-25 (Rowland et al. 2024) and CRISTAL-05 (Posses et al. 2025). For REBELS-25, we estimated Σ_SFR and Σ_[CII] assuming a 2D Gaussian distribution and the effective radius from Rowland et al. (2024) for estimating the area. For the A1689-zD1 components, we also assumed a 2D Gaussian distribution with the model FWHM from Table 3. We note that the plotted errors do not include the uncertainties of the sizes. Given the large uncertainties on the SFR estimates, the individual components of A1689-zD1 are comparable to the high-redshift results of REBELS-25, CRISTAL-05, and the relation from ALPINE. The red component appears partly consistent with the local relation from Herrera-Camus et al. (2015), which could be indicative of a more efficient heating (e.g., Herrera-Camus et al. 2015).

Fig. 11. Surface brightness, SFR, [C II], ΣSFR, and Σ[CII] for the individual components of A1689-zD1 compared to recent results of the resolved high-z galaxies REBELS-25 (Rowland et al. 2024) and CRISTAL-05 (Posses et al. 2025, the circles show resolved measurements across the galaxy). The turquoise line and region shows the local relation from Herrera-Camus et al. (2015), the blue dotted line and region show the relation from ALPINE (Schaerer et al. 2020, assuming an average area of 8 kpc2). The dashed gray line and region show the relation for z = 5.5 − 7 galaxies (Carniani et al. 2018).

4.1.3. SFR–L_[OIII]/L_[CII] relation

We also compared the line luminosity ratio of [C II] and [O III] to the SFR for z > 6 galaxies, shown in Fig. 12. While the estimates of the line ratios of the components span a range of a factor two, they are relatively consistent with the overall result for A1689-zD1. It has been found that L_[OIII]/L_[CII] is typically above 2 for z > 6 galaxies, however, this is based on a relatively small number of sources so far (Harikane et al. 2020; Carniani et al. 2020). High L_[OIII]/L_[CII] ratios have also been found in local low-metallicity dwarf galaxies (e.g., Cormier et al. 2015). However for local luminous IR galaxies, the ratio tends to be ≤1 (e.g., Howell et al. 2010; Díaz-Santos et al. 2017).

Fig. 12. Line luminosity ratio of [O III] and [C II] vs. SFR for each component of A1689-zD1 compared to other high-z detections (Inoue et al. 2016; Laporte et al. 2017, 2019; Carniani et al. 2017; Hashimoto et al. 2018, 2019; Tamura et al. 2019; Harikane et al. 2020; Witstok et al. 2022; Algera et al. 2024; Bakx et al. 2024). The diamonds represent the individual components color-coded according to their name, and the star formation is the one derived from the aperture measurements, as given in Tables 6 and 5. The orange and gray region and line indicate the range of the line ratio for local low-metallicity dwarf galaxies and starburst galaxies from De Looze et al. (2014).

Using a compilation of new data and literature results, Harikane et al. (2020) modeled the line ratios. They use simple photoionization models and discuss how several different physical effects impact the ratio; for example, the metallicity, ionization parameter, gas density, C/O ratio, and cosmic microwave background (CMB) attenuation. They find different properties can explain this, in particular very high ionization parameters or low photo-dissociation region coverage. Extensive modeling using cosmological radiation-hydrodynamics simulation have been published by Katz et al. (2022). Those simulations show that several effects, such as ionization parameter, CMB attenuation, and the ISM gas density, can play a significant role;however, in order to reproduce the observed results, the C/O abundance need to be lower than the solar value (e.g., Bakx et al. 2024).

Comparing to the line ratios of the SHINING survey for local star-forming galaxies (Herrera-Camus et al. 2018), we note that the [O III]-to-[C II] luminosity ratio is about a factor two higher. Given the relatively well-sampled global spectral energy distribution (SED; Bakx et al. 2021), we estimate the continuum ratio S_63 μm/S_122 μm ∼ 1.2 − 1.7. In the SHINING sample, the sources with S_63 μm/S_122 μm > 1 are typically also the sources with L_[OIII]/L_[CII] > 1; however, we note that SHINING only has few sources with such high line ratios.

4.1.4. Size-velocity relation

In Fig. 13 we also compare the relation between size and velocity dispersion of the individual components with other spatially-resolved measurements from local to high-redshift galaxies. For size, we adopted the effective radius for the comparison, assuming for the Gaussian profiles that the half-width-at-half-maximum is comparable with the effective radius. By assuming virialized clouds in the same manner as Spilker et al. (2022), we also calculated the expected size and velocity dispersion at given mass surface density (illustrated by the dotted lines in Fig. 13).

Fig. 13. Radius-velocity relation of the clumps in A1689-zD1 compared to giant molecular clouds (GMCs) in the Milky Way (Heyer et al. 2009), in the nearby spiral M64 Rosolowsky & Blitz (2005), a compilation of other low redshift galaxy GMCs (Bolatto et al. 2008), in the so-called Cosmic Snake at z = 1.0 (Dessauges-Zavadsky et al. 2019), and in the clumpy, lensed SPT 0311−58 at z = 6.9 (Spilker et al. 2022). Lines of constant surface mass density are shown, following Spilker et al. (2022). The clumps in A1689-zD1 have a high surface mass density, 1.5 dex higher than low-z GMCs and similar to SPT0311−58 and the Cosmic Snake, though with clump sizes that are factors of several smaller than in the more massive galaxy SPT 0311−58 and much larger than observed in less star-forming galaxies at lower redshift. Our smallest components are unresolved, indicating the possible presence of even smaller clumps, possibly corresponding in size to those in the Cosmic Snake or M64.

We find that the individual components are different from the distribution of the clumps observed in the local galaxies, but fall between the parameter spaces observed in the high-redshift, lensed, star-forming galaxies of the “Cosmic Snake” at z = 1.036 (Dessauges-Zavadsky et al. 2019) and SPT 0311−58 at z = 6.900 (Spilker et al. 2022). The distribution of the individual components are aligned with the dotted line, showing that their mass surface densities are comparable with those of the star-forming clumps observed in the other high-redshift lensed galaxies. This might indicate that the internal star-forming mechanisms resemble each other, even though there is a large difference in the total SFRs of these high-redshift lensed galaxies, from a few tens to thousands of solar masses per year.

4.1.5. Kinematics

To analyze the kinematics of each component, we created a fiducial 3D data cube for each component by subtracting the 3D models of the other components (Sect. 3.5) and produced moment-0 and moment-1 maps. In Fig. 14, we show the moment-0 and moment-1 maps of both the [C II] and [O III] lines for each component. The moment-0 maps are shown in terms of their S/N, with the rms indicated in the figure caption. Two versions of the moment-1 are displayed: one showing the velocity relative to the overall galaxy (z = 7.133) and one with zero velocity centered in the component.

We find in the [C II] moment-1 map that blue component shows a symmetric velocity structure, from blueshifted in NE to redshifted in SW around its intensity peak. Apart from the blue component, we do not find clear velocity gradients in the moment-1 maps of any other component in the [C II] or [O III]. We thus regarded blue as a potential rotator and the other components as dispersion-dominated systems. We note, however, that blue does not show a similar velocity gradient in the [O III] moment-1 map. As discussed in Rizzo et al. (2021), this may be explained by the difference in the motions between ionized and cold neutral gases, where the ionized gas is more likely tracing turbulent motions affected by stellar feedback for example, although the rotation has been seen in [O III] at z ∼ 9 (Tokuoka et al. 2022). In fact, Akins et al. (2022) report different line profiles between [C II] and [O III] at the highest SFR surface density region, mostly corresponding to green, which may be explained by the outflow of the ionized gas.

4.1.6. Dynamical masses

We also estimated the dynamical masses, M_dyn, of the individual components. For blue, we assumed a disk-like gas distribution, given the potential rotation (Sect. 4.1.5), and followed Equation (4) in Dessauges-Zavadsky et al. (2020):

$$\begin{array}{c}\hfill \left(\frac{M_{\mathrm{dyn}}}{M_{\odot }}\right)=1.16\times {10}^5{\left(\frac{v_{\mathrm{circ}}}{\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2\left(\frac{2r_{\mathrm{e}}}{\mathrm{kpc}}\right),\end{array}$$(1)

where r_e is the effective radius, v_circ is the circular velocity of the gaseous disk, including the inclination correction, given by v_circ = 1.763σ₀/sin(i) and i = cos⁻¹(FWHM₂/FWHM₁).

For the other components, we assume dispersion-dominated systems and virialized bodies, following Equation (10) of Bothwell et al. (2013), as

$$\begin{array}{c}\hfill \left(\frac{M_{\mathrm{dyn}}}{M_{\odot }}\right)=1.56\times {10}^6{\left(\frac{{\sigma }_0}{\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2\left(\frac{r_{\mathrm{e}}}{\mathrm{kpc}}\right),\end{array}$$(2)

where σ is the line-of-sight velocity dispersion and r_e is the effect radius of the sources as described from the fitting above, and thus giving the dynamical mass within the effective radius. Because the [C II] and [O III] lines could trace different gas motions and regions, we estimated two M_dyn values for each component separately – [C II] (M_dyn, [CII]) and [O III] (M_dyn, [OIII). We note that given the complicated morphology and kinematics of the system, each component might not be virialized. Nonetheless, the estimates provide a meaningful estimate of the dynamical mass. We summarize the various estimates of M_dyn in Table 7.

4.2. Mass of A1689-zD1

While A1689-zD1 was originally discussed as a modest, star-forming galaxy (Bradley et al. 2008; Watson et al. 2015), there are now several lines of evidence that indicate the galaxy is much larger and more massive than the stellar and gas mass derived from the continuum emission. The 12 kpc extent of the [C II] emission (regardless of its origin; Akins et al. 2022) as well as the increased IR SFR (Bakx et al. 2021; Akins et al. 2022, and Sect. 2) point to a bigger system. Furthermore, the dynamical mass estimates of the clumps (Sect. 7) indicate that the galaxy may be much more massive than originally believed, with a summed dynamical mass about 1.5 × 10¹⁰ M_⊙. This is an order of magnitude larger than the estimated stellar mass based on SED modeling including Spitzer-IRAC 3.6 μm and 4.5 μm detections. This implies an unfeasible high gas mass fraction given the approximately solar metallicity inferred from the FIR emission line ratios (Killi et al. 2023). The total metal mass of the galaxy would then exceed what could be produced the core collapse supernovae in a Chabrier initial mass function. On this basis, the gas mass fraction of this large dynamical mass is unlikely to be more than about 50%. There is certainly sufficient dust in the galaxy: a screen of half the estimated dust mass in A1689-zD1 distributed uniformly across 2 kpc² would result in an A_V ∼ 40 mag using the Milky Way dust extinction to column density ratio (Watson 2011). Since the star-formation obscuration fraction is about 90% (Akins et al. 2022), the stellar mass is unlikely to be underestimated by more than this, namely, by more than an order of magnitude. The stellar mass is therefore likely to be about 1–2 × 10¹⁰ M_⊙. These issues are discussed in greater depth in Killi et al. (2023).

4.3. Considering the physical origin of A1689-zD1

The complex gas structure and kinematics of A1689-zD1 revealed by the ALMA observations are different and exhibit a higher degree of irregularity compared to other z > 7 galaxies, where the multiplicity is often characterized by two or three [C II] components (e.g., Matthee et al. 2017; Carniani et al. 2017; Hashimoto et al. 2019), with separations on scales of 2–7 kpc. Other observations show either [C II] that is separated from other components (e.g., [O III]) and extended (e.g., Maiolino et al. 2015; Marrone et al. 2018). Single-components at scales of about 3 kpc with some degree of rotation have also been claimed in UV-bright galaxies at z > 6.5 (e.g., Smit et al. 2018; Tokuoka et al. 2022). Other [C II] detections either do not have sufficient S/N or angular resolution to resolve the emission (e.g., Knudsen et al. 2016; Bradač et al. 2017; Pentericci et al. 2016).

Many recent results for galaxies at z > 4, particularly from ALMA observations, show very clear velocity gradients (e.g., Jones et al. 2017; Smit et al. 2018; Neeleman et al. 2020; Tokuoka et al. 2022; Rowland et al. 2024; Posses et al. 2025). In many cases, the observed velocity gradients are interpreted as gas rotation and, in some cases, this can also be aptly described by rotational models (e.g., Carniani et al. 2013; Litke et al. 2019; Jones et al. 2021; Rizzo et al. 2021). These results are based both on non-lensed sources and gravitationally lensed galaxies and also across a range of different luminosities and masses. Certainly, the same system can change appearance due simply to inclination angle or evolution over a relatively short time. Only a few hundred million years separating z = 7 from z = 5 − 6 would be sufficient for a system to change from disk rotation to merger-dominated, according to the simulations of Kohandel et al. (2019).

It is worth noting that our angular resolution is about 5–20 times better than galaxies studied by Smit et al. (2018) and Matthee et al. (2017) due to the observational setup and the gravitational lensing magnification, which allows us to clearly rule out a simple ordered rotation for A1689-zD1. While a somewhat lower resolution would likely be sufficient to demonstrate that A1689-zD1 is not fully explained based on simple rotation, the sensitivity of data we have here, combined with a factor ∼9.6 in gravitational magnification reveal structures that are rarely detected in galaxies at high redshifts. Our conclusion here is somewhat at odds with Wong et al. (2022), who argued for a clear velocity SE-NW across the galaxy, with the NW being more dust-obscured. However, the presence of the violet component co-spatial with the redder part of the galaxy complicates a simple merger picture, as does the possible rotation and extent of the blue component. The results for the lensed galaxy SPT 0311−58 at z = 6.9 observed with ALMA, has a comparable S/N and sub-kiloparsec resolution and also displays a clumpy structure (Spilker et al. 2022). This suggests that a high S/N combined with spatial resolution significantly below the kpc scale may be important in fully exploring the dynamical state of at least some z > 6 galaxies.

4.3.1. An early galaxy assembly possibly consisting of five small clumps

One possible interpretation of the detection of several components within projected distances of ∼3 kpc and the detection of extended [C II] emission is that we are witnessing the early assembly of the inner kiloparsec-scale region of a galaxy. In such a scenario each of the components would represent ISM/mass substructure, possibly as (super)giant molecular clouds or young stellar clusters in the assembly of the bulge region. Given the large relative offsets and velocities, the ordered motion seen in disk galaxies would not have been established thus far. The extended [C II] detection (Akins et al. 2022), would either be halo gas or possibly part of extended ISM, part of a bound system. If this was the case, we might expect an extended stellar component following the structure of the extended [C II], which could be detectable with JWST. As argued above, the stellar mass may be much larger than indicated by the Spitzer observations. However, simulations of early galaxies suggest that they are unlikely to have effective radii larger than 2–3 kpc (Roper et al. 2022), indicating that such an extended galaxy interpretation is less favored. This suggests that if A1689-zD1 has a large obscured stellar mass, it is most likely concentrated in the inner few kiloparsecs.

4.3.2. Galaxy merger interpretation

The blue, green, and red components dominate the galaxy’s emission. As discussed above (and essentially indicated by Wong et al. 2022), it is possible that the blue and green components could be associated with the two rest-frame UV components detected with HST. Given the velocity difference of 80 km s⁻¹ and projected physical separation of 0.8–1.4 kpc between these two (the range is given by the slightly different positions in the [C II] and [O III] data), one interpretation could be interaction and merger. The red component shows an extended tail-like structure, which could suggest that this is an additional component being pulled towards the green and blue, while the elongated structure represents a tidal tail of the gas. The red component is separated by 0.7–1 kpc in projection and 220 km s⁻¹ from the blue component and by 1.8–2.1 kpc in projection and 140 km s⁻¹ from the green. This could suggest a scenario of a merger event between green and blue, with red being a gas-rich structure or a scenario of a triple-merger, where the red component contains a more dust-obscured stellar population and a tidal disruption of the gas distribution. The sizes and distribution could be consistent with the ‘merger’ phase observed in the simulations of Kohandel et al. (2019). However, without the spatial information from the radial direction, such merger scenarios are difficult to constrain.

4.3.3. Combined scenario

In the two subsections above, we describe the assembly of the inner few kiloparsecs of a galaxy or a merger scenario of two or three galaxies. A third plausible interpretation could be a combination of the two. Given the complex, almost chaotic structure, blue, green, and yellow appear to align in one direction and could represent the inner region of a galaxy, while red and violet could comprise more galaxy components that are interacting and possibly being accreted onto the main galaxy structure. The much higher dust obscuration of these NW components may support such a scenario.

Even though the ALMA data is of exceptional quality and reveals substructure so far unmatched in high-z galaxies, two pieces of information are missing: the radial spatial distribution, which cannot be measured, and the distribution of the stellar mass. This latter may well be obtained through deep observations with JWST. Zoom-in numerical simulations are necessary to provide additional constraints, similar to those of Kohandel et al. (2019).

4.4. Different star-formation indicators in the galaxy

With an ionization potential of 35 eV, double ionized oxygen is typically seen in only in HII regions of young massive stars (or near AGNs), so the [O III] 88 μm line traces ionized gas (H II regions) and is vastly less dust extinguished than optical or UV tracers (e.g., [O III] 5007 Å). Thus, it should serve as an excellent tracer of recent star-formation. The first ionization potential of carbon is only 11.26 eV, below atomic hydrogen and well below that of double ionized oxygen, but just above the lowest Lyman-Werner band energy of 11.15 eV required to two-step photo-dissociate molecular hydrogen. Furthermore, C⁺ is thus both found in the ionized and neutral gas phases. Recent work using absorption spectroscopy suggests that [C II] can be used to trace the neutral atomic gas component when metallicity is accounted for (Madden et al. 2020; Heintz et al. 2021). The lines of [C II] and [O III] can thus be used to trace different parts of the gas associated with the star-forming regions and the diffuse ISM. At z > 6, a few detections of both lines have been obtained (e.g., Inoue et al. 2016; Laporte et al. 2017, 2019; Carniani et al. 2017; Hashimoto et al. 2018, 2019; Tamura et al. 2019; Harikane et al. 2020), revealing a large scatter in the ratios. While it is not possible with these data alone to estimate the density and metallicity of the gas, with sufficient additional data, such estimates may be made from ISM models (e.g., Killi et al. 2023).

To compare the projected UV emission with the FIR emission in detail, we plot in Fig. 15 ratio maps of the [O III], [C II] lines, as well as the continua near these lines, and the combined UV emission from the sum of the F125W, F140W, and F160W HST filter images. This gives us nine ratio maps showing the resolved, normalized line ratios across the galaxy.

Fig. 15. Ratios of the galaxy’s UV, [C II]158 μm, [O III]88 μm, and dust continuum emission. The images are convolved to the smallest common beam: the [O III] beam where those data are involved and the [C II] beam otherwise. Parts of the ratio maps with S/N below 3 in the denominator are masked white, with increasing transparency from 3 to 6. The color range is normalized across the S/N > 6 part of the map for the panels involving the HST data. Contours for the [C II] and [O III] line emission and the UV emission are shown as dotted black, dashed green, and solid violet lines.

Three star-forming tracers are present: the [O III] line, the UV continuum, and the dust continuum emission. The UV continuum is susceptible to dust extinction. However, the tracers also have different overall lifetimes, with the [O III] emission being the most immediate, requiring 35 eV photons to ionize the gas, while the UV emission is longer lived, and the dust heating even longer. The most notable feature of these ratio maps is that the UV emission comes only from the eastern side of the galaxy, while the dust and FIR line emission is much more extended, revealing an entire half of the galaxy to the west that has no detectable UV emission. Secondly, the [C II] emission also clearly comes from a much more extended area than the other FIR emission, a point explored in more detail in Akins et al. (2022). As noted above, the dusty, FIR-dominated side of the galaxy seems to correspond to the red and yellow velocity components in the emission line data. The yellow component is particularly unusual in that it is dust-continuum dominated; namely, the dust continuum dominates over all other emission components, including [C II]. This component, however, is mixed in 88 μm and 158 μm continuum emission, indicating that it is not clearly colder or warmer, for example, than the dust in the rest of the galaxy. It does, however, also have relatively weak [O III] emission. At first glance, this could suggest that the oldest star-forming population in the galaxy does in fact reside here.

While we have argued that the two major components observed in the HST data seem to correspond primarily to the blue and green components identified from the ALMA line emission, it is not a precise match and suggests a greater complexity than what can be resolved with the current data (see Fig. 7). The central green component emits strongly in all four tracers, suggesting a young star-forming knot, while the blue component appears to be more UV dominated, suggesting that this component is somewhat older and is least dust-rich at its northern-most extent.

Finally, the comparison with the high-resolution data of Spilker et al. (2022) and Dessauges-Zavadsky et al. (2019) show that A1689-zD1 has a very high mass surface density. We have argued that the majority of the inferred mass in these surface densities must be stellar mass rather than gas mass, so these results are not quite comparable. However, we expect that a significant fraction of the mass is made up of gas. The high surface density and the range of clump sizes, from unresolved to nearly kiloparsec scale, indicate that A1689-zD1 has some aspects in common with the most massive star-forming galaxies at high-z. None of our data so far show any indication of AGN activity in the galaxy (as opposed to the possible interpretation from⁴Wong et al. 2022).

5. Conclusions

We present GBT band Q and ALMA band 3, 6, and 8 observations of the gravitationally lensed, epoch of reionization galaxy, A1689-zD1. The deep ALMA observations combined with the large gravitational magnification, provides some of the deepest high angular resolution data for a z > 7 galaxy to date, enabling us to resolve the ISM kinematics down to scales of 200 pc.

• Two emission lines have been detected and identified as [C II] 158 μm and [O III] 88 μm, consistent with a redshift of z = 7.1332 ± 0.0005. The spectroscopic redshift is offset by about 13 000 ± 7000 km s⁻¹ from the continuum break identified in previous X-shooter data.

• The line luminosities are L_[CII] = (5.8 ± 0.1)×10⁸ L_⊙ and L_[OIII] = (1.3 ± 0.02)×10⁹ L_⊙ after correcting for gravitational magnification. The line luminosities relative to the SFR are consistent with Local Universe galaxies, suggesting the ISM excitation conditions are normal.

• The galaxy system is complex with at least five individual components. Based on iterative modeling we estimate the structure and kinematical properties of each component. The components appear to be interacting, hinting at a galaxy in its early assembly. The complex kinematics and morphology do not show signs of a single rotating disk. The complexity is in agreement with predictions of early galaxy assembly from recent numerical zoom-in simulations.

• Several components of the galaxy are entirely obscured in the rest-frame UV and might also be similarly obscured in the optical. This could significantly affect the stellar mass estimate, indicating that A1689-zD1 might be significantly more massive than previously believed. Deep rest-frame NIR observations are the most direct route to resolving thisquestion.

Overall, A1689-zD1 still appears to be something of a mystery, with kinematics and gas surface densities that indicate it is in the process of assembly, but the dust and stellar properties argue for a much older system, with major discrepancies between its dynamical and stellar masses and between its Lyα break and FIR line redshifts, along with a small, compact core, especially in the UV emission; in addition, there is the vast [C II], and to a lesser extent, dust halo to consider. JWST is certain to resolve at least some of these issues and discrepancies.

The resolved [C II] and [O III] kinematics imply a more complex diversity among z > 6 galaxies. High-angular-resolution observations of other similar galaxies are needed to improve our understanding of early assembly and fully characterize the early evolutionary stages of normal galaxies.

Acknowledgments

We thank the anonymous referee for their helpful comments and suggestions. We acknowledge support from the Nordic ALMA Regional Centre (ARC) node based at Onsala Space Observatory. The Nordic ARC node is funded through Swedish Research Council grant No 2019-00208. We thank the Nordic ARC node staff for very helpful discussions. KK acknowledges support from the Swedish Research Council (2015-05580) and the Knut and Alice Wallenberg Foundation (KAW 2017.0292). DW, SF, and NB are supported in part by Independent Research Fund Denmark grant DFF-7014-00017. The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant number 140. TB acknowledges support from the Knut and Alice Wallenberg Foundation (KAW 2020.0081). A.K.I. is supported by NAOJ ALMA Scientific Research grant No.2020-16B. MJM acknowledges the support of the National Science Centre, Poland through the SONATA BIS grant 2018/30/E/ST9/00208. This research was funded in part by the National Science Centre, Poland, grant number 2023/49/B/ST9/00066. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.01406.S and 2017.1.00775.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

References

All Tables

All Figures

Fig. 1. Velocity-integrated maps of [C II] (left) and [O III] (right) of A1689-zD1 with the moment-0 (top) and moment-1 (bottom), derived from the velocity ranges v[CII] = −215 − 270 km s−1 and v[OIII] = −250 − 270 km s−1, respectively. The contours represent the 3, 5, 8, 12, and 17σ levels, where σ[CII] = 11.0 mJy km s−1 beam−1 and σ[OIII] = 63.2 mJy km s−1 beam−1. The moment-1 map is masked at 3σ of the moment-0 map.

In the text

Fig. 2. Mosaic of the [C II] binned by velocity in 40 km s−1 intervals. The value in the upper right corner of each sub-image indicates the middle of the velocity relative to z = 7.133. The gray contours show the band 6 continuum, where the levels correspond to 4, 6, 8σ. The five components of the analysis (Sect. 3.5) are indicated by name. The color bar below the last panel shows the level corresponding to the velocity integrated flux, SdV, for each map.

In the text

Fig. 3. Mosaic of the [O III] observations. The upper right corner shows the central velocity of each 40 km s−1 channel relative to z = 7.133. The gray contours show the band 8 continuum, where the levels correspond to 4, 6, 8, 10, 12, and 14σ. The five components of the analysis (Sect. 3.5) are indicated by name. The color bar below the last panel shows the level corresponding to the velocity integrated flux, SdV, for each map.

In the text

	Fig. 4. Non-detection of CO(3–2) emission measured at the GBT smoothed to a channel resolution of about 30 km s−1. The expected line location based on the [C II] and [O III] redshift of z = 7.133 is shown by the vertical dotted line.
In the text

	Fig. 5. [C II] (top) and [O III] (bottom) spectra (continuum-subtracted) overlaid with the resulting models for each component (color coded according to their component names). The black line shows the sum of the models. The residuals after subtracting the model results are shown as a thin gray histogram.
In the text

Fig. 6. Overall moment-0 maps (first and third rows) and moment-1 maps (second and fourth rows) for A1689-zD1 for [C II] and [O III] for the data (left), the model (middle), and the residuals (right). Contours show the model with the gray showing the moment-0 (3, 4, 5, 7, 10, 13, 16, and 19σ, where the rms is measured in the data regions excluding the source), and the white contours showing the moment-1 in steps of 20 km s−1. The applied mask is determined from the model. We note that this does not include the extended “halo” component presented in Akins et al. (2022).

In the text

	Fig. 7. Contours of the individual components of the [C II] 158 μm (top) and [O III] 88 μm (bottom). Overlaid is the combined HST image obtained with F125W, F140W, and F160W filters. The contours are color-coded according to their component name; contours are in steps of 2σ starting at 4σ (see Fig. 14 for the rms).
In the text

	Fig. 8. Decomposed components of the observations of A1689-zD1 shown in 3D with color-coding according to component name. The integrated moment-0 maps is shown at the end of the plot.
In the text

Fig. 9. Decomposition of the band 6 (about 158 μm, top) and band 8 (about 88 μm, bottom) dust continuum for the five components. The 4″ × 4″ cutouts of the observed, IMFIT model, and residual maps are presented from left to right. The solid and dashed contours show the intensity at the 2, 3, 4, 5, 6, 8, and 10σ levels and −3 and −2σ levels, respectively, which are measured based on the residual map. The spatial position of the five [O III] components are marked with the color crosses. The cyan circles indicate the elliptical apertures that are used for the aperture photometry.

In the text

Fig. 10. Line luminosity vs. SFR for L[CII] (left) and L[OIII] (right). Black outline diamonds with color centers indicate the color-coded individual components from Sect. 3.5. The uncertainty on the SFR-axis for each component corresponds to the range given by a Tdust = 25 − 55 K. The orange and gray lines and associated regions show the relations for star-forming galaxies and low-metallicity dwarf galaxies from De Looze et al. (2014). Data from the literature is included for comparison for high-z galaxies without AGNs (Capak et al. 2015; Inoue et al. 2016; Knudsen et al. 2016; Pentericci et al. 2016; Bradač et al. 2017; Laporte et al. 2017, 2019, 2021; Carniani et al. 2017; Smit et al. 2018; Walter et al. 2018; Hashimoto et al. 2018, 2019; Tamura et al. 2019; Harikane et al. 2020; Binggeli et al. 2021; Fujimoto et al. 2021; Molyneux et al. 2022). We note that the data from the REBELS survey include the sources with L[CII] and UV + IR SED-estimated SFR (Sommovigo et al. 2022; Ferrara et al. 2022). On the right panel, for [O III] we also include the results from Marrone et al. (2018). In the left panel, the blue and magenta lines represent model predictions from Vallini et al. (2015) and Olsen et al. (2017).

In the text

Fig. 14. Moment-0 and moment-1 maps for each component: [C II] results (left) and [O III] results (right). The models of the other components have been subtracted from the data cube. The left columns show the moment-0 maps integrated over the line of the component; the middle columns the moment-1 given relative to the A1689-zD1 overall velocity field with zero velocity corresponding to z = 7.133; and the right columns the moment-1 map for each individual velocity (same velocity used for both [C II] and [O III]). The overlaid contours represent the fit models, in the moment-0 maps these correspond to 3, 4, 5, 6, 7, 8, 9, and 10σ; in the moment-1 maps, these correspond to velocity in km s−1 in steps of 10 km s−1. The moment-0 maps are given in S/N, where the rms in each map is measured excluding the source region (for [C II]: rms = 5.4, 11.6, 11.5, 6.4, and 8.3 mJy km s−1 beam−1, respectively, and for [O III]: rms = 27.7, 50.8, 54.1, 29.4, and 38.1 mJy km s−1 beam−1, respectively (from top to bottom).

In the text

Fig. 11. Surface brightness, SFR, [C II], ΣSFR, and Σ[CII] for the individual components of A1689-zD1 compared to recent results of the resolved high-z galaxies REBELS-25 (Rowland et al. 2024) and CRISTAL-05 (Posses et al. 2025, the circles show resolved measurements across the galaxy). The turquoise line and region shows the local relation from Herrera-Camus et al. (2015), the blue dotted line and region show the relation from ALPINE (Schaerer et al. 2020, assuming an average area of 8 kpc2). The dashed gray line and region show the relation for z = 5.5 − 7 galaxies (Carniani et al. 2018).

In the text

Fig. 12. Line luminosity ratio of [O III] and [C II] vs. SFR for each component of A1689-zD1 compared to other high-z detections (Inoue et al. 2016; Laporte et al. 2017, 2019; Carniani et al. 2017; Hashimoto et al. 2018, 2019; Tamura et al. 2019; Harikane et al. 2020; Witstok et al. 2022; Algera et al. 2024; Bakx et al. 2024). The diamonds represent the individual components color-coded according to their name, and the star formation is the one derived from the aperture measurements, as given in Tables 6 and 5. The orange and gray region and line indicate the range of the line ratio for local low-metallicity dwarf galaxies and starburst galaxies from De Looze et al. (2014).

In the text

Fig. 13. Radius-velocity relation of the clumps in A1689-zD1 compared to giant molecular clouds (GMCs) in the Milky Way (Heyer et al. 2009), in the nearby spiral M64 Rosolowsky & Blitz (2005), a compilation of other low redshift galaxy GMCs (Bolatto et al. 2008), in the so-called Cosmic Snake at z = 1.0 (Dessauges-Zavadsky et al. 2019), and in the clumpy, lensed SPT 0311−58 at z = 6.9 (Spilker et al. 2022). Lines of constant surface mass density are shown, following Spilker et al. (2022). The clumps in A1689-zD1 have a high surface mass density, 1.5 dex higher than low-z GMCs and similar to SPT0311−58 and the Cosmic Snake, though with clump sizes that are factors of several smaller than in the more massive galaxy SPT 0311−58 and much larger than observed in less star-forming galaxies at lower redshift. Our smallest components are unresolved, indicating the possible presence of even smaller clumps, possibly corresponding in size to those in the Cosmic Snake or M64.

In the text

Fig. 15. Ratios of the galaxy’s UV, [C II]158 μm, [O III]88 μm, and dust continuum emission. The images are convolved to the smallest common beam: the [O III] beam where those data are involved and the [C II] beam otherwise. Parts of the ratio maps with S/N below 3 in the denominator are masked white, with increasing transparency from 3 to 6. The color range is normalized across the S/N > 6 part of the map for the panels involving the HST data. Contours for the [C II] and [O III] line emission and the UV emission are shown as dotted black, dashed green, and solid violet lines.

In the text