Ionised gas kinematics and dynamical masses of z ≳ 6 galaxies from JADES/NIRSpec high-resolution spectroscopy

We explore the kinematic gas properties of six 5 . 5 < z < 7 . 4 galaxies in the JWST Advanced Deep Extragalactic Survey (JADES), using high-resolution JWST / NIRSpec multi-object spectroscopy of the rest-frame optical emission lines [O iii ] and H α . The objects are small and of low stellar mass ( ∼ 1kpc; M ∗ ∼ 10 7 − 9 M ⊙ ), less massive than any galaxy studied kinematically at z > 1 thus far. The cold gas masses implied by the observed star formation rates are ∼ 10 × larger than the stellar masses. We find that their ionised gas is spatially resolved by JWST, with evidence for broadened lines and spatial velocity gradients. Using a simple thin-disc model, we fit these data with a novel forward modelling software that accounts for the complex geometry, point spread function, and pixellation of the NIRSpec instrument. We find the sample to include both rotation-and dispersion-dominated structures, as we detect velocity gradients of v ( r e ) ≈ 100 − 150kms − 1 , and find velocity dispersions of σ 0 ≈ 30 − 70kms − 1 that are comparable to those at cosmic noon. The dynamical masses implied by these models ( M dyn ∼ 10 9 − 10 M ⊙ ) are larger than the stellar masses by up to a factor 40, and larger than the total baryonic mass (gas + stars) by a factor of ∼ 3. Qualitatively, this result is robust even if the observed velocity gradients reflect ongoing mergers rather than rotating discs. Unless the observed emission line kinematics is dominated by outflows, this implies that the centres of these galaxies are dark-matter dominated or that star formation is 3 × less e ffi cient, leading to higher inferred gas masses.


Introduction
In the nearby Universe galaxies show a variety of dynamical structures and structural components, which are reflective of their mass assembly histories (e.g.Cappellari 2016;van de Sande et al. 2018;Falcón-Barroso et al. 2019).However, the details of the formation and evolution of these structures -nominally, rotationally supported discs and spheroidal bulges supported primarily by dispersion -are still unclear.The physical conditions in the early Universe, the secular evolution of galaxies, and mergers with other systems are all likely to play important roles.One outstanding question in particular is when and how early galaxies settled into dynamically cold discs.
Although this question is ideally answered by measuring spatially resolved stellar kinematics across cosmic time, such measurements have only been possible up to z ∼ 1 (e.g, van Houdt et al. 2021), except for a few strongly lensed, massive galaxies at z ∼ 2 (Newman et al. 2018).Instead, the ionised gas of the interstellar medium (ISM) provides critical insight into the dynamical properties of (star-forming) galaxies across a much wider redshift range (for a review, see Förster Schreiber & Wuyts 2020).Many studies have focused on inferring galaxy dynamical properties from rest-frame optical emission lines, to map the evolution in the velocity dispersion (σ) and the ratio between ⋆ degraaff@mpia.de the rotation velocity and dispersion (v/σ), which measures the degree of rotational support of the system.
Measurements of the ionised gas kinematics from multiple large spectroscopic surveys of star-forming galaxies at z ∼ 1 − 4 have demonstrated that the velocity dispersion of star-forming galaxies increases with redshift, while the rotational support decreases to v/σ ≈ 1 by z ≈ 3 (e.g., Wisnioski et al. 2015Wisnioski et al. , 2019;;Stott et al. 2016;Simons et al. 2017;Turner et al. 2017;Price et al. 2020).In addition, imaging studies have shown that galaxy morphologies are less disc-like and more clumpy at rest-frame UV wavelengths at higher redshifts (e.g., van der Wel et al. 2014;Guo et al. 2015;Zhang et al. 2019;Sattari et al. 2023).Theoretical models and simulations have suggested that gravitational instabilities, the accretion of gas and smaller systems from the cosmic web, and stellar feedback may be responsible for increased turbulence at higher redshifts (e.g., Dekel et al. 2009b;Genel et al. 2012;Ceverino et al. 2012;Krumholz et al. 2018).
However, using submillimeter observations from the Atacama Large Millimeter Array (ALMA), several studies have found dynamically cold discs at z ∼ 2 − 6 (e.g., Neeleman et al. 2020;Jones et al. 2021;Lelli et al. 2021;Rizzo et al. 2021;Parlanti et al. 2023;Pope et al. 2023), even finding v/σ ≈ 20 at z ≈ 4.5 (Fraternali et al. 2021).These observations raise the question of how such systems formed and settled so rapidly (within ∼ 1 Gyr), and how these observations can be reconciled with the aforementioned studies at cosmic noon.How-ever, the ALMA observations infer the galaxy kinematics from far-infrared and millimetre transitions (CO, [Cii]), which trace colder gas than the rest-frame optical lines, likely explaining part of the discrepancy (Übler et al. 2019;Rizzo et al. 2023).Additionally, selection effects likely play an important role, as many of the ALMA observations primarily probe the most massive galaxies at z > 4.
To understand the evolution of more typical (∼ M * ) galaxies, requires spatially resolved spectroscopy of faint galaxies at high redshifts.In this regime, ground-based telescopes are unable to observe rest-frame optical emission lines, whereas sub-mm facilities are in principle able to observe such systems, but at extremely high cost.The launch of the James Webb Space Telescope (JWST) has enabled spectroscopy with very high sensitivity and high spatial resolution (Gardner et al. 2023;Rigby et al. 2023).Using the slitless spectroscopy mode of JWST/NIRCam (Rieke et al. 2023b), Nelson et al. (in prep.)reveal the ionised gas kinematics in a massive galaxy at z ≈ 5.However, only the NIRSpec instrument provides the spectral resolution needed to resolve galaxy kinematics for low-mass systems (σ ≈ 50 km s −1 for a uniformly-illuminated slit; Jakobsen et al. 2022).JWST/NIRSpec additionally provides a multi-object spectroscopic (MOS) mode (Ferruit et al. 2022), allowing for the simultaneous observation of up to ≈ 200 objects, making observations of high-redshift targets highly efficient.The slit-based observations with the microshutter array (MSA) however sacrifice one spatial dimension with respect to integral field spectroscopy (IFS; Böker et al. 2022, see also Price et al. 2016).Therefore, extra care is required to extract spatial and dynamical information from NIRSpec MSA data.
In this paper, we present the dynamical properties of six high-redshift galaxies (z > 5.5) in the JWST Advanced Deep Extragalactic Survey (JADES; Eisenstein et al. 2023).These objects are spatially extended in deep NIRCam imaging and were followed up with the high-resolution NIRSpec MOS mode, providing spatially-resolved spectroscopy of their rest-frame optical emission lines.The data are presented in Section 2. To model the galaxy kinematics, we propagate analytical models through a simulated NIRSpec instrument, and use MCMC sampling to fit the data (Section 3), using NIRCam imaging as a prior on the morphology.We present the results of our modelling in Section 4, demonstrating a diverse range of kinematic structures in a previously unexplored population of galaxies.We discuss the possibility that some systems may be late-stage mergers in Section 5, and examine the large discrepancy between the derived dynamical masses and stellar masses.We summarise our findings in Section 6.

NIRSpec spectroscopy
We use NIRSpec MOS observations in the GOODS-South field taken as part of the JADES deep and medium programmes (ID numbers 1210 and 1286, PI N. Lützgendorf;Bunker et al. 2023;Eisenstein et al. 2023).Targets were selected from a combination of JWST/NIRCam and Hubble Space Telescope (HST) imaging and followed up with JWST/NIRSpec using the low-resolution prism (R ∼ 100), the three medium-resolution gratings (R ∼ 1000), and the reddest high-resolution grating (G395H; R ∼ 2700 for a uniformly illuminated slit).Here, we focus primarily on the high-resolution spectroscopy, although we also use the prism data to estimate stellar masses and star formation rates (SFRs).The spectra in our sample were obtained using a 3-point nodding pattern and vary in depth, ranging from 2.2 h to 7.0 h of total integration time for the G395H grating (summarised in Table 1).The total exposure times for the prism range from 2.2 h to 28 h.
The NIRSpec data were reduced using the NIRSpec GTO collaboration pipeline (Carniani et al. in prep), as is also described in Curtis-Lake et al. (2023) and Bunker et al. (2023).Crucially, in contrast with other studies so far using NIRSpec data, our analysis largely does not rely on the final rectified, combined, and extracted 2D and 1D spectra generated by the pipeline.Instead, we perform our dynamical modelling to intermediate data products: we use 2D cutouts of the detector from individual exposures that have been background-subtracted and flat-fielded.In this way, we mitigate correlated noise and artificial broadening that is otherwise introduced by the rectification and combination algorithm of the reduction pipeline.We also note that these intermediate data products do not correct for any slit losses due to the spatial extent of the sources, as these effects are already fully accounted for in our modelling (Section 3).
Nevertheless, the 2D rectified high-resolution spectra and their 1D extractions were used for the initial visual inspection and selection of the sample.We also used the 1D extracted prism data for spectral energy distribution (SED) modelling in Section 4 to estimate stellar masses and SFRs.

NIRCam imaging
s Even though some objects were initially selected based on HST imaging, JWST/NIRCam imaging is available for all targets in our sample from a combination of Cycle 1 programmes.The majority of our targets fall within the JADES footprint in GOODS-S (Rieke et al. 2023a), and are therefore imaged in 9 different NIR-Cam filters.One of the selected targets (JADES-NS-10016374) is located outside of the JADES footprint but falls within the footprint of the FRESCO survey (programme 1895, PI P. Oesch;Oesch et al. 2023).Although FRESCO is primarily a grism survey, the survey also obtained imaging with three different NIR-Cam filters (F182M, F210M, F444W), although at significantly reduced exposure times compared with JADES imaging.
All images were reduced as described in Rieke et al. (2023a).We use the JWST Calibration Pipeline v1.9.6 with the CRDS pipeline mapping (pmap) context 1084.We run Stage 1 and Stage 2 of the pipeline with the default parameters, but provided our own sky-flat for the flat-fielding.Following Stage 2, we perform a custom subtraction of the 1/f noise, scatteredlight effects ("wisps") and the large-scale background.We perform an astrometric alignment using a custom version of JWST TweakReg, aligning our images to the HST F814W and F160W mosaics in the GOODS-S field with astrometry tied to Gaia-EDR3 (G.Brammer priv.comm.).We achieve an overall good alignment with relative offsets between bands of less than 0.1 short-wavelength pixel (< 3 mas).We then run Stage 3 of the JWST pipeline, combining all exposures of a given filter and a given visit.
For our analysis we select the NIRCam filter that most closely represents the emission line morphology of the target.Given the high equivalent widths of the emission lines in our sample, we use the medium band covering the emission line where available (four out of six objects).For the remaining two objects we instead use broad-band filters.We use the available low-resolution prism spectra to quantify the flux originating from emission lines versus the stellar continuum in Appendix B. For the four objects with medium-band images, we hence estimate that ≈ 70 − 75% of the NIRCam flux is due to emission  .11572-27.77495 JADES-NS-10016374 53.11572 -27.77496 12 16.8 Hα JADES-GS+53.18374-27.79390 JADES-NS-20086025 53.18375 -27.79389 9 8.0 [Oiii] Notes. (a) The ID number corresponds to the NIRSpec ID as described in Bunker et al. (2023). (b) Best-fit coordinates from the morphological modelling to the NIRCam imaging.These differ slightly from the coordinates in the JADES ID due to updated astrometry.
lines, and therefore provide a good map of the ionised gas.For the two objects with broad-band images the stellar continuum dominates (emission line fluxes contribute ≈ 35%), and we discuss how this may affect the inferred kinematic parameters in Sections 4 and 5.2.

Sample selection
We have visually inspected all (358) JADES targets in GOODS-S for which high-resolution NIRSpec spectroscopy is available thus far.We select objects that are at high redshift (z > 5.5, i.e. when the Universe was less than 1 Gyr old), are spatially extended in the NIRCam imaging (r e ≳ 0.1 ′′ ), and have bright emission lines, i.e. an integrated signal-to-noise ratio (SNR) ≳ 20.Additionally, we require that the 1D spectrum shows no obvious evidence for a broad-line component, which would be indicative of large-scale outflows (i.e., the sample discussed in Carniani et al. 2023).The resulting sample consists of six objects that span the redshift range z = 5.5 − 7.4.We observe Hα in two of these objects, and [Oiii] in the remaining four.For the two highest redshift targets Hα is outside the wavelength coverage of NIRSpec (z > 7); for two z ∼ 6 objects we do not observe Hα, as the traces of the high-resolution spectra are long and therefore fall partially off the detector.The 2D combined and rectified images of the emission lines are presented in Fig. 1 together with cutouts from the NIRCam imaging, where we have selected the NIRCam filter closest to the emission line as described in the previous section.The positions of the microshutters are shown in orange.
We note that, given the complex selection function of JADES (Bunker et al. 2023), and our additional imposed selection criteria during the visual inspection, the selected sample is far from complete in (stellar) mass, magnitude or star formation rate.However, our aim is to demonstrate the ability of NIRSpec/MOS to measure galaxy kinematics in an entirely new regime: Fig. 2 shows that the targets in our sample are not only at a higher redshift than has been attainable so far with ground-based nearinfrared spectrographs but are also significantly smaller and less massive.We defer a comprehensive analysis of the galaxy kinematics as a function of redshift and mass to a future paper, as this will require the complete JADES and NIRSpec WIDE datasets (Eisenstein et al. 2023, Maseda et al. in prep.) as well as a thorough understanding of the selection functions of the different survey tiers.

Dynamical modelling
We use a Bayesian forward-modelling approach to estimate the dynamical properties of the galaxies from 2D emission spectra, observed through the NIRSpec MSA apertures.First, we construct parametric model cubes for the flux distribution (I(x, λ)) based on analytical surface brightness and velocity profiles.Second, we model the complexities of the NIRSpec instrument that are imprinted on the data when mapping the kinematic models onto a mock NIRSpec detector.Third, we use a Markov chain Monte Carlo (MCMC) sampling method to fit the models to the spectra, adopting Sérsic profile fits to NIRCam imaging as a prior.We defer a detailed description of this forward-modelling and fitting software (msafit) to a future paper, in which we will also demonstrate convergence tests and comparison with calibration data.In this section, we only provide a summary overview of the models and software, which we release publicly together with this paper1 .

Thin disc models
Although the spatial and spectral resolution of JWST/NIRSpec are high, the small systems in our sample are close to the resolution limit.This suggests that we should limit ourselves to a relatively simple geometric and dynamical model: a thin rotating disc.Although geometrically thin, we allow for the disc to be kinematically warm by adding a velocity dispersion profile.We discuss the possible limitations of our model choice in Section 5.
We model the spatial flux distribution of the emission line as a Sérsic profile (Sérsic 1968), which is described by four parameters: the total flux F, half-light radius r e (major axis), projected minor-to-major axis ratio q, and Sérsic index n.As we assume a thin disc model, the projected axis ratio is directly related to the inclination angle (i) of the system.In addition, there are three important position-dependent parameters that enter the model, the position angle (PA) with respect to the MSA slitlet as measured from the positive x-axis (i.e.PA = 90 deg represents perfect alignment with the slit), and the centroid position of the object within the shutter (dx, dy).
For the velocity field we use the common empirical description of an arctangent rotation curve (Courteau 1997): where v a is the asymptotic or maximum velocity with respect to the systemic velocity of the galaxy, and r t is the turnover radius.
We parametrise the systemic velocity as the mean wavelength of the emission line λ 0 .To allow for a kinematically warm disc, we additionally assume a constant velocity dispersion profile σ(r) ≡ σ 0 across the disc.
Combined, this amounts to a 11-parameter model (λ 0 , F, dx, dy, r e , n, q, PA, r t , v a , σ 0 ).We convolve the flux pro-file and velocity profiles to form a model flux cube I(x, y, λ), where x and y are the coordinates in the MSA plane and sampled in user-specified intervals.To construct these cubes, we ensure that the spatial and wavelength dimensions are sampled at minimum at Nyquist frequency for the point spread function (PSF) at the wavelength considered or the NIRSpec pixel size (0.1 ′′ ), whichever is smaller.However, this sampling would be too sparse to evaluate the steep Sérsic profiles at small radii.In order to accurately integrate Sérsic profiles we therefore first oversample the spatial grid dynamically, such that the innermost region (< 0.2 r e ) is oversampled by a factor 500 and the outer regions (> r e ) by a factor 10, and then integrate the profile onto the coarser grid.

Forward modelling the NIRSpec MSA
Although forward-modelling software for slit-based multi-object spectroscopy has been developed before (Price et al. 2016), there are several unique challenges to modelling NIRSpec MOS data: (i) the diffraction-limited PSF of JWST, which enables high spatial resolution, but is highly complex in shape; (ii) the complex geometry of the NIRSpec MSA, comprising ≈ 2 × 10 5 microshutters that are separated by shutter walls, imprints additional diffraction patterns (+-shaped due to the slit aperture) as well as shadows on the detector ("bar shadows"); (iii) the relatively large pixels (≈ 0.1 ′′ × 0.1 ′′ ) of the NIRSpec detector imply that the PSF is undersampled at all wavelengths 2 .
As a consequence of the first and second above challenges, there is a strong variation of the PSF shape within the shutter.High-redshift extragalactic objects are often also comparable in size to the NIRSpec PSF width and pixel size, and the position of the flux centroid within the shutter therefore strongly affects the shape of the flux distribution on the detector.Moreover, the shutter walls (≈ 0.035 ′′ ) are relatively small compared to the pixel size, making the effects of the bar shadows complex to model.Lastly, the open area of the microshutters is relatively small (≈ 0.2 ′′ × 0.46 ′′ ; Jakobsen et al. 2022) compared to the PSF size.Slit losses are therefore substantial, even at the centre of the shutter.
We attempt to capture all these effects within our modelling.First, we construct libraries of synthetic PSF models at a range of different wavelengths and spatial offsets, with the PSF centres sampling the shutter every 0.02 ′′ and using a 5x oversampling factor for the PSF images.These PSFs represent the 2D image of a point source with an infinitely narrow emission line in the detector plane, and hence contains both the spatial distribution along the slit and the distribution in the wavelength direction.We constructed the PSFs using custom Fourier optical simulations, tracing monochromatic point sources through NIRSpec to the detector focal plane.These models capture the combined PSF of JWST and NIRSpec, including the diffraction and light losses (often referred to as path losses) caused by the masking by the micro-shutter slitlets and spectrograph pupil.We defer a detailed description to de Graaff et al. (in prep.), but note that the construction of these PSFs is largely the same as in previous works that presented or used the NIRSpec Instrument Performance Simulator (Piquéras et al. 2008(Piquéras et al. , 2010;;Giardino et al. 2019;Jakobsen et al. 2022).The main difference is that the implementation used for this work is python-based and makes use of the Physical Optics Propagation in python (POPPY; Perrin et al. 2012) libraries, allowing a carefully tuned wavelength-independent sampling in both the image and pupil planes.Although these models are based on in-flight calibrations where possible, a number of necessary reference files were created prelaunch.We therefore caution that there is likely to be a systematic uncertainty in the true width and shape of these PSFs.Unfortunately, neither sufficient nor dedicated calibration data currently exist.We discuss the current status of calibrations in more detail in Appendix A, and estimate a ≈ 10 − 20% systematic uncertainty in the PSF full width at half maxima (FWHM; both spatial and spectral) of our models.
Second, we construct libraries of spectral traces for all microshutters and dispersers, using the instrument model of Dorner et al. (2016) and Giardino et al. (2016), the parameters of which were tuned during the in-flight commissioning phase (Lützgendorf et al. 2022;Alves de Oliveira et al. 2022).These traces provide a mapping from the centre of a given shutter (s i j ) and chosen wavelength to the detector plane (X, Y).We also use this model to derive the tilt angle of the slitlets with respect to the trace in the detector plane (a few degrees for G395H).
Third, we construct model detectors, for which the pixels are initially oversampled by a factor 5. To reduce computational cost we do not model the full 2048×2048 detectors, but create cutouts of ≈ 30 × 30 pixels around a region of interest, whilst keeping track of the corresponding detector coordinates.
With these libraries and models in place, we can forward model the analytical flux cube of Section 3.1.As the PSF strongly varies with the intrashutter position, the model cube cannot be convolved with a single PSF.Instead, we treat the model cube as a collection of point sources, hence propagating each point in the cube with its local PSF.The slices of the model cube are then projected onto the oversampled model detector using the trace library, given a defined shutter (s i j ).Finally, the detector is downsampled by a factor 5 to match the true detector pixel size, and convolved with a 3 × 3 kernel to mimic the effects of inter-pixel capacitive coupling3 .This results in a noiseless model image of the input data cube on the detector.

Model fitting
The procedure of Section 3.2 generates a model for a single set of parameters.To estimate the posterior probability distributions of the parameters, we use the MCMC ensemble sampler implemented in the emcee package (Foreman-Mackey et al. 2013).
Importantly, we perform the comparison between the models and data in the detector plane to mitigate correlated noise.We hence do not use the 2D combined spectrum, but simultaneously fit to multiple exposures (Fig. 3), while masking pixels that are flagged as being affected by cosmic rays or hot pixels.The likelihood function for a set of parameters p model then is where N is the number of exposures, K the number of unmasked pixels per exposure, and F j , M j and σ j are the observed flux, model flux, and uncertainty in the j th pixel, respectively.
In calculating the posterior, we use informative priors where possible, as the geometry is poorly constrained based on the  1).Although the final combination of all exposures (left) was used for our initial visual inspection and sample selection, the pixels in this spectrum are highly correlated.Instead of using this combined spectrum, we simultaneously fit to all individual exposures obtained.In the case of JADES-NS-00016745 two exposures were taken per nodding position, resulting in six independent measurements for one 3-point nodding pattern with NIRSpec.To combat the undersampled PSF of NIRSpec, we perform our modelling in the detector plane, propagating parametric models to the exact same location on the detector as the observed data.The likelihood is then computed from the combination of all residual images.Pixels flagged by the reduction pipeline as affected by cosmic rays are masked and shown in grey.
spectroscopic data alone.We perform morphological modelling to the NIRCam images (Fig. 1) with lenstronomy (Birrer & Amara 2018;Birrer et al. 2021), following the procedure described in Suess et al. (2023).Based on the SNR of the object in the image and corresponding typical uncertainties in the structural parameters derived by van der Wel et al. ( 2012), we set Gaussian priors centred on the best-fit lenstronomy estimates.To allow for uncertainty in the PSF models and deviations between the image morphology and emission line morphology (as described in Section 2.2; Appendix B), we double all uncertainties in the structural parameters to set the dispersions of the Gaussian priors.The mean wavelength and integrated flux of the emission line can be determined from the 1D spectrum with high accuracy.This line flux needs a correction for the slit losses, which we estimate based on the best-fit lenstronomy model parameters.We then create Gaussian priors for the central wavelength and flux, somewhat conservatively assuming an uncertainty of half a pixel (≈ 2Å) for the wavelength, and a 10% uncertainty in the total flux on the detector.Lastly, we allow for a small uncertainty on the intrashutter position of the source due to the finite pointing accuracy of JWST, for which we use Gaussian priors with a dispersion of 25 mas, which is the typical pointing accuracy after the MSA target acquisition (Böker et al. 2023).
For the parameters of the dynamical model, the maximum velocity (v a ) and the velocity dispersion (σ 0 ), we use uniform priors.For the fitting we parametrise the turnover radius as the ratio r t /r e , and assume a uniform prior for this ratio.We show an example fit in Figs. 3 and 4, which demonstrates that v a and σ 0 are formally well-constrained (> 5 σ significance).The turnover radius is poorly constrained, and forms the largest source of uncertainty on v a due to the degeneracy between v a and r t /r e .This is likely due to a combination of the moderate spatial resolution and the limited spatial extent probed by the microshutters.In Section 4 we instead compute the rotational velocity at r e , v(r e ), which is better constrained by the data (blue distributions in the top right of Fig. 4).
We note that the relatively large spatial extent of the source compared to the shutter size may also lead to a loss of flux, as the background subtraction step in the reduction pipeline subtracts flux from the (neighbouring) source that falls in the adjacent shutter.To test the magnitude of this potential bias, we also perform our modelling on a separate reduction that excludes exposures in which the source falls in the central shutter and includes only the outer nods, therefore mitigating any selfsubtraction and contamination (but at the cost of a slightly lower overall SNR).We find that the recovered flux is consistent within the error bars with the fit to the standard reduction that uses all 3 nodding positions.Our model is therefore robust against slight self-subtraction present in the spectra, likely helped by the prior information provided by the NIRCam imaging and the relatively small spatial extent of the sources compared to the shutter size.

Results
4.1.Ionised gas kinematics at z > 5.5 We present the results of our modelling for the sample of six objects in Tables 2 and 3. Significant rotation is detected in three of the objects, and marginally detected or consistent with zero in the other cases.As we did not remove galaxies that are strongly misaligned with the slit, some of these objects (e.g., object JADES-NS-00019606) may have a velocity gradient that is simply not observable at this position angle of the MSA.Nevertheless, the measurement of the velocity dispersion is still useful in these cases, and also provides confirmation that our model is able to return v a ≈ 0 despite our a priori assumption that the system is rotating.
We find five objects have velocity dispersions that are broader than the instrument line-spread function (LSF) for a point source (σ inst ≈ 25 − 30 km s −1 , see Appendix A), which is the relevant LSF after accounting for the source morphology.The formal uncertainties on the dispersions are small, which is likely due to our assumption of a thin disc, meaning that our estimate of the error on σ 0 does not include an uncertainty caused Extrapolating the fit by Übler et al. (2019) at z ∼ 1 − 3 to z ∼ 7 suggests that the ionised gas in galaxies is highly turbulent at early epochs.In contrast, we find that all objects lie well below this extrapolation, and instead have velocity dispersions that  are approximately equal to the average dispersion at z ∼ 2 − 3.
On the other hand, the typical stellar mass of our sample is substantially lower than the literature data (Fig. 2).If the velocity dispersion depends on stellar mass (as predicted by simulations, e.g.Pillepich et al. 2019), the ISM in these low-mass systems may still have a relatively high turbulence.
We also compare with studies that used ALMA to resolve galaxy kinematics at the same redshift as our sample (blue triangles; Neeleman et al. 2020;Rizzo et al. 2020;Fraternali et al. 2021;Lelli et al. 2021;Rizzo et al. 2021;Herrera-Camus et al. 2022;Parlanti et al. 2023).Although these objects lie at the same redshift, the measurements differ substantially: the galaxies observed with ALMA are often more massive (M * ≳ 10 10 M ⊙ ), and the observed emission lines tend to trace much colder gas.Interestingly, despite these differences, the ALMA-based velocity dispersions are very similar to our measurements of the ionised gas based on rest-frame optical emission lines.Possibly, this is because the effects of the higher mass and lower gas temperature on the velocity dispersion act in opposite directions.Observations of the same systems with both ALMA and JWST will be crucial to constrain these effects.
Next, we compute the ratio v/σ ≡ v(r e )/σ 0 and examine its dependence on redshift in Fig. 6, comparing with the same literature as mentioned previously.Studies around cosmic noon showed a clear, gradual decline in the degree of rotational sup- .Rotational support as a function of redshift, measured as the ratio between the velocity at the effective radius and the constant velocity dispersion: v(r e )/σ 0 .Although studies based on ground-based nearinfrared data (as described in Fig. 5) have found a clear, gradual decline in v/σ toward higher redshift, we find an interesting diversity among our sample of low-mass galaxies, with dynamically-cold discs existing possibly as early as z ≈ 7.
port toward higher redshift.Based on these measurements, one may expect none of the z > 5 galaxies to be rotation-dominated (v/σ > 1).Yet, we find an interesting diversity among our sample, with three objects having v/σ > 1 even at the highest redshifts (z ≈ 7).We discuss in Section 5 whether these objects may truly form cold rotating discs, or whether these reflect velocity gradients within systems that are not virialised.
We again compare our sample with the ALMA-based studies, which are all rotation-dominated systems with relatively high v/σ ratios.Our sample shows greater diversity, which may be due to the fact that the gas tracers differ and the mass range probed is significantly different.The misalignment of some objects with the microshutters also may underestimate the intrinsic v/σ ratio of some systems.Larger samples are therefore required to fully understand the different kinematic properties of the gas phases traced with ALMA and JWST at high redshifts Lastly, we revisit Section 2.2, where we described that for two out six objects the NIRCam image used as a prior in the emission line modelling predominantly traces stellar continuum emission instead of line emission.If the morphology of the emission line differs strongly from the continuum, this may bias the inferred kinematic parameters, especially if the galaxy is prolate instead of our assumed oblate thin disc model.One of these two objects (JADES-NS-00016745; Fig. 3) has a major axis in the NIRCam image that is well-aligned with the microshutter, and we observe a strong velocity gradient in the 2D spectrum with only a small offset between the major axis PA from the imaging and the median of the posterior distribution of the PA (shown in Fig. 4).A prolate morphology is therefore highly unlikely for this object.The nature of the second object (JADES-NS-100016374) is more uncertain, as we only marginally detect rotation.If the kinematic major axis of the ionised gas differs strongly from the photometric major axis, the true rotational velocity and v/σ 0 may be substantially higher than inferred with our modelling.

Comparing dynamical and stellar masses
We use the dynamical models to examine the mass budget of the galaxies.For a system in virial equilibrium, the dynamical mass enclosed within radius r is computed as where v circ is the circular velocity, k the virial coefficient, and G the gravitational constant.However, for comparison with the total stellar mass, we define a 'total' dynamical mass as described in Price et al. (2022): As we assume a thin disc model, we adopt k tot = 1.8, which is the virial coefficient for an oblate potential with q = 0.2 and n ∼ 1 − 4 (Price et al. 2022, Fig. 4).The true shape of the potential is not well-constrained however, and this choice for k tot therefore introduces a systematic uncertainty in the dynamical mass estimates, as k tot can vary by up to a factor two.Following Burkert et al. (2010), we compute the circular velocity as which accounts for the effects of pressure gradients on the rotational velocity, and depends on the disk scale length (r d ).
At the effective radius the pressure correction term reduces to 2(r e /r d ) = 3.35.We note that for the one object with uncertain oblate/prolate morphology (Section 4.1), this calculation of v circ may be incorrect.However, as discussed in Section 5.2, the inferred dynamical mass is likely less affected.
Next, we compare these total dynamical masses to stellar masses.To estimate stellar masses and SFRs, we perform SED modelling with the Bayesian fitting code BEAGLE (Chevallard & Charlot 2016) to the low-resolution prism spectra.The fits were run adopting a two-component star formation history consisting of a delayed exponential with current burst, a Chabrier (2003) initial mass function with an upper mass limit of 100 M ⊙ , and a Charlot & Fall (2000) dust attenuation law assuming 40% of the dust in the diffuse ISM.We note that the 1D prism spectra were flux-calibrated assuming a point-like morphology and without considering NIRCam photometry.Although this slit loss correction approximately corrects for the variation in the PSF FWHM with wavelength, there is a systematic offset between the total flux of the object and the flux captured by the slit.We estimate this aperture correction using our modelling software and the morphology in the long wavelength filter (F444W; measured with lenstronomy as described in Section 3.3), finding correction factors in the range 1.2 − 2.5, and apply this to the stellar masses and SFRs.The inferred properties are presented in Appendix B, together with an example prism spectrum and SED model.
We compare the estimated dynamical and stellar masses in Fig. 7, and for reference plot the same ground-based nearinfrared studies as in Figs. 5 and 6.As may be expected from the fact that M dyn includes dark and baryonic mass, all objects in M d y n = M * Fig. 7. Stellar mass versus dynamical mass (Eq.4) as inferred from the prism and high-resolution spectroscopy, respectively.The dashed line shows the one-to-one relation between the two masses.Data points from the literature (circles, squares) are as described in Fig. 5.As is to be expected, the dynamical masses are larger than the stellar masses for all objects in our sample.Surprisingly, however, the dynamical masses are substantially larger (up to a factor ≈ 40), most likely indicating large gas masses or large systematic uncertainties in the stellar mass estimates.
our sample have dynamical masses that are greater than the estimated stellar masses.However, the difference between the two masses is much larger than in previous studies, on average deviating by as much as a factor 30.Only Topping et al. (2022) have reported similarly large stellar-to-dynamical mass discrepancies at z ∼ 7, albeit for more massive galaxies and based on spatiallyintegrated line width measurements instead of spatially-resolved dynamical modelling as in this paper.We discuss the possible origins of this discrepancy between the stellar and dynamical masses in detail in Section 5.2.

Discussion
These data and modelling have taken us to a very new regime of galaxy kinematics: low-mass galaxies at z > 5. Our modelling results are formally very well constrained, and the resulting parameter constraints (e.g., M * vs. M dyn ) at face value imply spectacular results.Yet, a look at Fig. 1 also makes it clear that our simple symmetric models may not capture the complex geometry of the systems.Therefore, our results warrant and require careful discussion.

Clumpy cold discs or mergers?
Using a forward modelling approach we have been able to separately constrain the morphology, velocity gradient and intrinsic velocity dispersion for each JADES object.To do so, we assumed an underlying model of a thin rotating disc (Section 3.1).
The velocity gradients measured in the context of our model suggest that the systems are dynamically relatively cold with higher than anticipated v/σ ratios (Fig. 6) based on an extrapolation of kinematic studies at z < 4.However, both observations and theoretical models have suggested that the rate of (major) mergers rises rapidly toward z ∼ 6  Notes. (a) The intrashutter offsets are measured in terms of shutter pitch rather than arcsec, as this unit is constant for all shutters and does not depend on the spatial distortion across the MSA. (b) The morphology of this source is also discussed in Baker et al. (2023).
Table 3. Dynamical modelling results: dynamical properties and model-derived quantities.Values are the median of the posterior probability distributions, and uncertainties reflect the 16 th and 84 th percentiles.Notes. (a) The sign of the model parameter v a indicates the observed direction of the velocity gradient along the slit.
(e.g., Rodriguez-Gomez et al. 2015;Bowler et al. 2017;Duncan et al. 2019;O'Leary et al. 2021).It is therefore likely that some objects in our sample are merging systems, or have recently merged with another galaxy.We indeed find complex (emission line) morphologies for several objects, most notably apparent in objects JADES-NS-00016745 and JADES-NS-20086025 (see Fig. 1), and Baker et al. (2023) show that object JADES-NS-00047100 can be described by three separate morphological components.It is therefore possible that the rotational velocities inferred under the assumption of a virialised system may actually reflect the velocity offset between two (or more) objects, or velocity gradients induced by the gravitational interaction in a pre-or post-merger phase.
Observations, on the other hand, also show that high-redshift galaxies often contain large star-forming clumps, and that the overall clumpiness of galaxies increases toward higher redshifts and lower masses (e.g., Guo et al. 2015;Carniani et al. 2018;Zhang et al. 2019;Sattari et al. 2023).Importantly, these clumps do not necessarily lead to a globally unstable system, and can be sustained within a rotationally-supported, albeit warm, disc (Förster Schreiber et al. 2011;Mandelker et al. 2014).
With the small angular scales (∼ 0.2 ′′ ) and velocity differences (∆v ∼ 200 km s −1 ) involved for the systems studied in this paper, it is very difficult to distinguish between a merging system and star-forming clumps with ordered rotation.This degeneracy has been discussed extensively in the literature, although for lower redshifts and larger angular scales (e.g., Krajnović et al. 2006;Shapiro et al. 2008;Wisnioski et al. 2015;Rodrigues et al. 2017).Simons et al. (2019) used simulations of merging galaxies to construct mock observations and hence quantify the frequency with which these systems are misclassified as rotating discs, showing that misclassifications are very common (≈ 50%), unless very stringent disc selection criteria are applied.Similarly, Hung et al. (2015Hung et al. ( , 2016) ) demonstrated that it becomes increasingly difficult to distinguish mergers from rotating systems toward later stages in the interaction between galaxies.On the other hand, Robertson et al. (2006) used hydrodynamical simulations to show that mergers between gas-rich systems can also lead to the formation of rotating discs with high angular momentum.Gravitational interaction between galaxies and the formation of rotating discs are thus not necessarily mutually exclusive, and the high gas fractions inferred in the next Section (5.2) are at face value consistent with the scenario proposed in these simulations.
Therefore, for any individual galaxy in our sample we cannot definitively conclude whether it is truly a rotating disc or an ongoing merging system.Although the NIRSpec MOS data are unprecedented in depth, resolution and sensitivity for galaxies at this mass and redshift, the objects are resolved by only a few resolution elements along a single spatial direction.Follow-up observations with the NIRSpec IFS mode can provide resolved 2D velocity field maps for these systems, which may then be compared against the disc selection criteria of Wisnioski et al. (2015) and Simons et al. (2019) to improve the constraints on their dynamical states.However, the high-resolution NIRSpec IFS observations needed are not feasible for large samples of objects.It therefore seems inevitable at present to accept the fact that the nature of individual galaxies remains ambiguous.A sta-tistical framework combining merger rates of galaxies and their observed emission line kinematics may provide a way forward to observationally constrain the settling of galaxies into cold discs at high redshifts, with number statistics that will be provided by upcoming surveys.

Uncertainties in the mass budget
We found a large discrepancy between the stellar and dynamical masses for the JADES objects.The dynamical mass presumably reflects the sum of the dark, stellar and gas mass within the effective radius.It is therefore not unexpected that the dynamical masses are larger than the stellar masses.However, the magnitude of the mass discrepancy -more than an order of magnitude -is surprising, as it is significantly larger than previous studies at lower redshifts at low mass (e.g., Maseda et al. 2013).
Considering the discussion of the previous section, we should first examine the robustness of our dynamical masses.To calculate the dynamical mass (Eq.4) we assumed that the galaxies are virialised; to compute the circular velocity, we assumed that the mass profile is approximately consistent with a rotating disc and exponential mass distribution.For the dispersiondominated objects, the latter assumption may be problematic.If we instead assume a spherical mass distribution for these objects, we can instead follow the dynamical mass calculation for dispersion-supported systems by Cappellari et al. (2006): where the virial coefficient depends on the Sérsic index, K(n) = 8.87 − 0.831n + 0.0241n 2 .However, for the low Sérsic indices of our sample K(n ≈ 1) ≈ 8 and is therefore comparable to the coefficients that enter Eq. 4, as k tot v 2 circ ≈ 3.35k tot σ 2 0 ≈ 7σ 2 0 .Similarly, in case of a prolate mass distribution we would expect k tot ∼ 4 (Price et al. 2022) and v circ ≈ √ 3σ 0 .In other words, the dynamical mass would not be overestimated by much if systems were actually dispersion dominated.On the contrary, for the object in our sample with uncertain oblate/prolate morphology (Section 4.1) the rotational velocity is possibly underestimated, which would lead to an underestimation of the dynamical mass and stellar-to-dynamical mass discrepancy.
On the other hand, for the rotation-dominated objects M dyn will be dominated by the value of v(r e ).If this velocity does not reflect a rotational velocity of a virialised system, but a velocity offset between two objects, then we cannot expect the dynamical mass estimate to be accurate.Both Simons et al. (2019) and Kohandel et al. (2019) explored the effects of incorrect physical and observational assumptions on the resulting v circ and M dyn estimates using mock observations of simulated galaxies.Simons et al. (2019) showed that for a merging system (noting that this is only a single simulation), the circular velocity is on average overestimated by a factor ≈ 1.5 (≈ 0.15 dex), which translates into an overestimation of M dyn by a factor 2 (0.3 dex).Kohandel et al. (2019) showed that, depending on the assumed inclination angle, M dyn can be both under-and overestimated in the case of a merger, and find a mass discrepancy of ≈ ±0.3 dex for velocity offsets of the same magnitude as found in our sample.Together with the uncertainty in the virial coefficient (see Section 4.2), we therefore conclude that the dynamical masses may be overestimated by up to ≈ 0.3 − 0.6 dex, which cannot explain the large differences we find between the stellar and dynamical masses.
Yet, in the above we assume that the inferred gas kinematics are dominated by gravitational motions.If the velocity . The baryonic, stellar and gas masses (estimated using the relation between Σ SFR and Σ gas from Kennicutt 1998) as a fraction of the dynamical mass.We find that the gas mass, and hence the baryonic mass, is approximately one order of magnitude larger than the stellar mass.Although the inclusion of the gas component reduces the large discrepancy in mass found in Fig. 7, a factor 3 − 4 difference between the dynamical and baryonic mass still remains for all but one object.
dispersions or velocity gradients are instead the result of nongravitational motions, i.e. turbulence and outflows due to stellar feedback, then the dynamical masses may be severely overestimated.As is discussed in great detail in Übler et al. (2019), based on theoretical models the turbulence due to stellar feedback appears to be in the range of ∼ 10 − 20 km s −1 (Ostriker & Shetty 2011;Shetty & Ostriker 2012;Krumholz et al. 2018).This is significantly lower than the circular velocities measured for our sample and therefore cannot lead to a large bias in our dynamical masses.Outflows, however, may form a larger source of uncertainty.We have selected against objects in JADES with outflows as presented in Carniani et al. (2023), who measured outflow velocities v out > 200 km s −1 .Lower outflow velocities are difficult to detect, but may still be present in our data.We therefore turn to observations of starburst galaxies at low redshift for comparison.Heckman et al. (2015) and Xu et al. (2022) detected outflows using UV metal absorption lines, and demonstrated that the ratio v out /v circ correlates with both the specific SFR (sSFR ≡ SFR/M * ) and the SFR surface density (Σ SFR ).Based on Fig. 10 of Xu et al. (2022) and the fact that for our sample sSFR ∼ 10 −8 yr −1 and Σ SFR ∼ 10 M ⊙ yr −1 kpc −2 , we estimate v out /v circ ∼ 3.This would imply an overestimation of the circular velocity by a factor 3 or a factor 10 in the dynamical mass.However, it is unclear whether the outflowing gas traced by the rest-frame UV absorption lines is also traced by the restframe optical emission lines, and how this in turn would translate into the uncertainty on the dynamical mass.For example, Erb et al. (2006) found no correlation between the Hα line width (and hence the dynamical mass) and galactic outflow velocities measured from rest-frame UV absorption lines for more massive star-forming galaxies at z ≈ 2.
If the dynamical mass is robust (at the factor 2 level), we should turn our attention to the other mass components that contribute to the total mass budget.An obvious component not dis-cussed so far is the gas mass.Both observational and theoretical studies have shown that the gas content is important to take into account (Price et al. 2016;Wuyts et al. 2016), as the typical gas fraction rises rapidly toward higher redshift and lower masses (for a review, see Tacconi et al. 2020).We estimate the gas masses of our sample based on the SFRs obtained from the SED modelling to the prism spectroscopy.We use the inverse of the Kennicutt (1998) relation between gas surface mass density (Σ gas ) and SFR to infer the total gas mass.Although calibrated only at low redshifts, this likely provides a reasonable order-ofmagnitude estimation for the high SFR surface densities of our sample (Daddi et al. 2010;Kennicutt & Evans 2012).
From this, we find gas masses of the order M gas ∼ 10 9 M ⊙ with the average ratio ⟨M gas /M * ⟩ ≈ 10.We can now compare these gas masses and resulting baryonic masses (M bar ≡ M * + M gas ) to the dynamical masses in Fig. 8.Although the gas mass is large compared to the stellar mass ( f gas ≡ M gas /M bar ≈ 0.90; consistent with measurements of more massive galaxies at z ∼ 7 from Heintz et al. 2022), a discrepancy of approximately a factor 3 − 4 between the baryonic and dynamical mass remains for all but one object.Our estimate of the gas masses carries a systematic uncertainty, as it for instance assumes a constant star formation efficiency.Price et al. (2020) show that a different estimator for the gas mass (following Tacconi et al. 2018) results in a mean gas mass difference of only 0.13 dex for a sample of more massive galaxies at cosmic noon.The same gas scaling relations are poorly, if at all, constrained at z ∼ 6 in the stellar mass regime considered in this work, and we therefore cannot make the same comparison for our sample.To increase the gas masses and make the baryonic masses consistent with the dynamical masses would, however, require a factor ≈ 3 decrease in the star formation efficiency for the majority of our sample.Yet, Pillepich et al. (2019) show that in the cosmological simulation TNG50 the gas fraction in low-mass galaxies (M * ∼ 10 9 ) increases rapidly with redshift up to z ∼ 3−4, but then appears to flatten at higher redshifts with M gas /M dyn ≈ 0.2.Further observations as well as simulated data at even lower stellar masses and higher redshifts will be necessary to constrain the gas masses for objects similar to the JADES sample.
Importantly, also the stellar masses may suffer from systematic uncertainties.The low-resolution spectroscopy spans a broad wavelength range (1 − 5 µm), probing the rest-frame UV to optical for all objects, and therefore provides good constraints on the recent star formation history (SFH).However, these measurements may suffer from an 'outshining effect' in which a young star-forming population dominates the SED, making it near-impossible to detect an underlying population of older stars at rest-frame optical wavelengths (e.g., Maraston et al. 2010;Pforr et al. 2012;Sorba & Sawicki 2018;Giménez-Arteaga et al. 2023;Tacchella et al. 2023;Whitler et al. 2023).This is especially relevant for our sample, as we selected bright emission lines to perform our dynamical modelling, and these lines tend to have high equivalent widths.Using mock observations of cosmological simulations, Narayanan et al. (2023) show that the outshining effect may underestimate the stellar masses by 0.1−1.0dex at z ≈ 7 depending on the selected prior for the SFH.This makes the possible severe underestimate of the stellar mass the most problematic potential source of systematic errors in our mass budget.Imaging at longer wavelengths with JWST/MIRI and spatially-resolved SED fitting may offer improvements in the stellar mass estimates in the future (e.g., Abdurro'uf et al. 2023;Giménez-Arteaga et al. 2023;Pérez-González et al. 2023).We note, however, that even this effect would play little role in the regime where M dyn is 30 × M * .
In summary, a number of systematic effects may be contribute to the large discrepancy of a factor 30 between the stellar and dynamical masses for this low-mass, high-redshift sample.We argue that the dynamical masses are relatively wellconstrained, with at most a 0.3 − 0.6 dex uncertainty even in the case of an ongoing merger, although we cannot rule out a bias in some of the dynamical mass measurements due to galactic-scale outflows.Clearly, a substantial amount of gas must be present in these highly star-forming galaxies, and we estimate this can account for ≈ 1 dex in the stellar-to-dynamical mass discrepancy.Although the gas masses are uncertain, and scaling relations between SFR and gas surface densities have significant scatter, we believe it is unlikely that the gas masses are underestimated by a large factor for all objects in our sample.This leaves the possibility that instead the stellar masses may be significantly underestimated, as we lack constraints at longer rest-frame wavelengths where an older stellar population may be measurable.
In this discussion we have neglected one mass component thus far: the dark matter.On the small spatial scales probed (≈ 1 kpc) this may not appear to be a dominant factor in the mass budget, particularly as multiple studies have shown a rapid increase in the central baryon fraction toward higher redshifts (e.g., van Dokkum et al. 2015;Price et al. 2016;Wuyts et al. 2016;Genzel et al. 2017Genzel et al. , 2020).Yet, it is interesting to consider a situation where there is significant dark matter within the effective radii at these redshifts.Within the ΛCDM cosmological model, dark matter dominates the mass content of the Universe.Under hierarchical structure formation, small dark matter haloes form first, whereas stars only form after the gas within those haloes has cooled sufficiently, subsequently growing into galaxies through accretion (cold gas streams, mergers; e.g., White & Rees 1978;Dekel et al. 2009a;Oser et al. 2010).In this scenario, it may be possible that very young galaxies are dark matter dominated even in the central regions, as the baryonic mass is still under assembly.This is of particular interest, as galaxies in the stellar mass regime probed in this paper (∼ 10 8 M ⊙ ) may be the progenitors of galaxies of M * ∼ 10 11 M ⊙ at z = 0 (Moster et al. 2018;Behroozi et al. 2019), which have baryon-dominated centres.As the aforementioned observations at cosmic noon are typically also of more massive galaxies (M * ∼ 10 10−11 M ⊙ ), the difference in the central baryon fractions between our work and measurements at cosmic noon may therefore not be in tension, but reflect the time evolution in the distribution of dark and baryonic mass as galaxies grow and assemble their stellar mass.We explore this idea further using cosmological simulations in de Graaff et al. (in prep.).However, to definitively conclude whether this scenario may apply to the JADES galaxies will require a thorough understanding of the systematic uncertainties on the different mass components.

Conclusions
We use the JADES spectroscopic sample in GOODS-S to select six targets at z = 5.5 − 7.4 that are spatially extended in NIRCam imaging, and for which high-resolution (R ∼ 2700) spectroscopy was obtained with the NIRSpec MSA.We show that these galaxies lie in a previously unprobed part of parameter space: not only because of their high redshifts, but also their small sizes (∼ 1 kpc) and low stellar masses (M * ∼ 10 8 M ⊙ ).The high-resolution spectra reveal rest-frame optical emission lines ([Oiii] and Hα) that are broadened and have spatial velocity gradients.
To extract the dynamical properties we model the objects as thin, but warm, rotating discs.We describe a novel forward mod-de Graaff et al.: Resolving galaxy kinematics at z > 6 with JWST/NIRSpec elling software to account for several complexities of data taken with the NIRSpec instrument: the PSF, shutter geometry and bar shadows, and pixellation.Using NIRCam imaging as a prior on the emission line morphology, we are able to constrain the rotational velocities and velocity dispersions of the objects in our sample, and hence also estimate dynamical masses.Our findings can be summarised as follows.
-The objects in our sample are small (r e ∼ 0.5 − 2 kpc), of low stellar mass (M * ∼ 10 7.5−8.9 ) and modest star formation rates (SFR ∼ 2 − 20 M ⊙ yr −1 ), which we infer from SED modelling to low-resolution NIRSpec spectroscopy.The gas masses implied by the SFRs are on average 10× larger than the stellar masses.-We find intrinsic velocity dispersions in the range σ 0 ≈ 30 − 70 km s −1 , which is consistent with studies reporting the velocity dispersions of more massive galaxies at cosmic noon.
-Three out of six objects show significant spatial velocity gradients, resulting in v/σ ≈ 1 − 6.Under the assumption of our thin disc model, this implies that the high-redshift objects are rotation-dominated discs.However, we cannot rule out the possibility that the detected velocity gradients reflect velocity offsets between interacting galaxies.-Comparison between the dynamical and stellar masses reveals a surprising discrepancy of a factor 10 − 40.After accounting for the large gas masses, the dynamical masses still remain larger than the baryonic masses by a factor ∼ 3. -We argue that the dynamical masses are robust within a factor 2 − 4 even in the case of an ongoing merger.Only the presence of outflows, if these were to dominate the observed emission line kinematics, can substantially lower the inferred dynamical masses.However, the baryonic-todynamical mass discrepancy might also imply that the centres of these objects are dark-matter dominated.Moreover, there are large systematic uncertainties on the stellar and gas masses.The baryonic masses can be reconciled with the dynamical masses if the star formation efficiency in these objects is a factor 3 lower than initially assumed.
Our work provides a first demonstration of the powerful capabilities of the NIRSpec MOS mode to perform spatially-and spectrally-resolved analyses.Crucially, this enables the study of galaxy kinematics in a highly efficient manner, as a single observation can probe a wide range in redshift, mass and SFR.With larger spectroscopic samples using the high-resolution MOS mode currently being acquired, JWST NIRSpec will in the nearfuture allow for statistical analyses of the origins and settling of disc galaxies in the early Universe.

Fig. 1 .Fig. 2 .
Fig. 1.Sample of six spatially-resolved high-redshift objects in JADES.Left panels show cutouts of the emission lines in the 2D rectified and combined spectra obtained with the high-resolution G395H grating.Negatives in the cutouts are the result of the background subtraction method used.Right panels show NIRCam image cutouts for each object (JADES, FRESCO), for the band that most closely resembles the emission line morphology (Section 2.2).The 3-shutter slits and 3-point nodding pattern used result in an effective area of 5 shutters: the shutter encompassing the source is shown in orange, and the shutters used for background subtraction are shown in purple.

Fig. 3 .
Fig.3.Example of the fitting procedure for object JADES-NS-00016745 (Fig.1).Although the final combination of all exposures (left) was used for our initial visual inspection and sample selection, the pixels in this spectrum are highly correlated.Instead of using this combined spectrum, we simultaneously fit to all individual exposures obtained.In the case of JADES-NS-00016745 two exposures were taken per nodding position, resulting in six independent measurements for one 3-point nodding pattern with NIRSpec.To combat the undersampled PSF of NIRSpec, we perform our modelling in the detector plane, propagating parametric models to the exact same location on the detector as the observed data.The likelihood is then computed from the combination of all residual images.Pixels flagged by the reduction pipeline as affected by cosmic rays are masked and shown in grey.

deFig. 4 .
Fig.4.Corner plot for the 11-parameter thin-disc model for the object in Fig.3.Histograms show the posterior probability distributions, with orange lines indicating the prior probability distributions.We generally find good constraints on the kinematic parameters v a and σ 0 , although the turnover radius r t is poorly constrained and degenerate with the rotational velocity.The top right panels show (in blue) the parameters that we derive from the model and are discussed in Section 4.
Fig.6.Rotational support as a function of redshift, measured as the ratio between the velocity at the effective radius and the constant velocity dispersion: v(r e )/σ 0 .Although studies based on ground-based nearinfrared data (as described in Fig.5) have found a clear, gradual decline in v/σ toward higher redshift, we find an interesting diversity among our sample of low-mass galaxies, with dynamically-cold discs existing possibly as early as z ≈ 7.

Table 1 .
Coordinates and G395H exposure times of the selected sample.

Table 2 .
Dynamical modelling results: morphological properties and wavelength.Values are the median of the posterior probability distributions, and uncertainties reflect the 16 th and 84 th percentiles.