The ALPINE-CRISTAL-JWST Survey: The Tully-Fisher Relation at $z = 4-6$

Introduction

Understanding how galaxies assembled their mass and structure in the early Universe is a key goal of extragalactic astrophysics. During the first Gyr of cosmic history, typical star-forming galaxies were rapidly growing through gas accretion and mergers, while their internal structure was still settling [63, 13, 56, 94, 43]. This early phase is key for studying how baryons assembled within dark matter (DM) haloes, and when rotationally supported galaxies began to emerge.

Gas kinematics provides a straightforward way to address the dynamical state of galaxies. For example, integral field unit (IFU) surveys observing the ionized-gas kinematics of galaxies at $z \approx1 - 3$ revealed a diverse galaxy population, including turbulent rotating disks, dispersion-dominated systems, and mergers [31, 117]. These galaxies were found, on average, exhibiting larger velocity dispersions and lower rotational support compared to their local counterparts [117, 36, 81]. However, the magnitude of this evolution remains debated and appears to depend on the gas tracer. Measurements based on cold gas generally show lower velocity dispersions and higher rotational support than ionized-gas measurements, implying a milder evolution of disk turbulence and rotational support with redshift [91, 95, 90] (Girard et al. 2021; Oliveira et al. 2023).

Extending such analysis to earlier epochs is crucial to determine when stable rotating disks become common, how disk settling proceeds with cosmic time, and how this process is connected to the build-up of stellar mass ($M_\star$) within galaxies.

In this context, a useful way to study the connection between the mass of a galaxy and its dynamics is the Tully-Fisher relation (TFR). Originally defined as a relationship between the luminosity and rotation velocity of systems [110], it is now commonly studied in terms of their stellar mass and/or baryonic mass ($M_{\rm bar}$) versus a characteristic circular velocity (e.g., [66, 88, 105, 61]). The stellar-mass TFR (sTFR) and the baryonic TFR (bTFR) provide a powerful framework for testing how disk-like galaxies assemble their mass within dark matter halos. The cosmic evolution of both the sTFR and bTFR can place constraints on gas fractions, star-formation efficiency, feedback, and the relative contributions of baryons and dark matter on galaxy scales (e.g., [112, 107, 5, 100]).

Observationally, the cosmic evolution of the TFR has been studied extensively up to $z \approx2$, with mixed results reported in the literature depending on the redshift range, sample selection, and galaxy population. Some studies report little or no evolution in the sTFR out to $z \approx1 - 2$, suggesting that the sTFR was already in place for massive star-forming disks at relatively early times [10, 67, 17, 107, 5]. In contrast, other studies report negative offsets with respect to the local relation, indicating that galaxies at fixed stellar mass show higher circular velocities at earlier cosmic times, possibly due to an evolution in the efficiency of gas-to-star conversion within growing dark matter halos [84, 115, 112, 41]. Overall, there is currently no clear consensus on the cosmic evolution of the sTFR. The disagreement among the studies can arise from observational and methodological differences in the derivation of galaxy kinematics, stellar masses, and data quality [112, 77, 100, 23, 64]. Compared to the sTFR, the bTFR cosmic evolution remains less constrained because robust measurements of gas masses are required, which remain challenging at high redshift. In the local Universe, adding the gas mass into the mass budget for dwarf galaxies moves these systems into a single baryonic mass TFR [66, 61]. At high redshift, accounting for gas mass becomes increasingly important for massive systems because galaxy gas fractions increase toward earlier cosmic times [106].

Probes of both the bTFR and the sTFR at $z \gtrsim4$ remain scarce, mainly due to observational limitations. Recent Atacama Large Millimeter/sub-millimeter Array (ALMA) observations of the [C II] 158 µm line emission have revealed that some galaxies at $z \approx4 - 6$ are consistent with being rotating disks [56, 49, 42, 58]. However, observations remain limited because long exposure times are required to properly spatially resolve the systems. The James Webb Space Telescope (JWST) has opened a new window on this problem. Early Near Infrared Camera (NIRCam) imaging studies showed that galaxies at $z \approx4 - 6$ span a wide range of morphologies, including disturbed systems, mergers, and a significant fraction of disk-like objects [28, 50, 57]. Some of the more massive systems even show bulges in the early stages of formation [45]. JWST IFU observations now make it possible to map rest-frame optical emission lines, such as Hα, on kpc scales at these epochs, allowing for a direct comparison between ionized-gas kinematics, stellar morphology, and cold-gas tracers. In parallel, similar multi-wavelength approaches are being extended even deeper into the end of the epoch of reionization [1] (Fisher et al. 2026). This new observational landscape makes it possible to explore the TFR at $z \gtrsim4$, in an epoch where galaxies are expected to be more gas-rich, more turbulent, and still assembling a significant fraction of their stellar and baryonic mass. Studying the TFR in this regime is therefore important because it probes the epoch when the connection between galaxy mass and gravitational potential may still be emerging, offering direct constraints on the early build-up of disks and on the relative role of baryons and dark matter in young galaxies.

In this paper, we use JWST/Near Infrared Spectrograph (NIRSpec) IFU observations towards a sample of main-sequence star-forming galaxies at $z = 4$–6, the ALPINE-CRISTAL-JWST survey [24]. The IFU data enable the characterization of the ionized-gas kinematics of galaxies at high resolution (0\farcs15 \sim1 kpc at $z = 5$), allowing more robust kinematic measurements than single-slit or slitless spectroscopy. Combined with JWST/NIRCam images, the JWST/NIRSpec IFU data allow us to identify disk-like systems. Our main goal is to study the Tully-Fisher relation in its baryonic and stellar-mass versions, and to explore possible cosmic evolution of the TFRs through comparisons with lower-$z$ samples (e.g., [88, 61, 112]).

The paper is organized as follows: Section 2 describes the sample, observations, and ancillary data used. In Section 3, we introduce the models adopted to analyze the data and our galaxy disk classification scheme. We present the baryonic and stellar-mass TFRs in Section 4 and discuss them in Section 5 by means of comparisons with other observational works and with a simple empirical model. Finally, we conclude in Section 6. We adopt a standard flat Λ-CDM cosmology with $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_{\mathrm{M}} = 0.3$, and $\Omega_{\Lambda} = 0.7$.

Observations and Data Reduction

We analyze JWST/NIRSpec IFU observations (Cycle-2 program GO–3045; P.I.: A. Faisst) of 18 main-sequence star-forming galaxies at $z = 4$–6, the ALPINE-CRISTAL-JWST survey. The targets were selected from the ALMA large programs ALPINE [56, 3, 26, 27] and CRISTAL [43], which provide [C II] observations at 0\farcs2–0\farcs7 spatial resolution ($\sim1$–5 kpc at $z = 5$). The ALPINE-CRISTAL-JWST survey benefits from the exquisite multi-wavelength imaging coverage of the COSMOS field [97, 116, 18, 101], which enables co-spatial comparisons with the stellar continuum observed at sub-kpc scales. This makes the ALPINE-CRISTAL-JWST sample a key reference to connect galaxy properties between the post-reionization era, cosmic noon, and the local universe. The sample spans stellar masses of $M_\star\approx10^{9-11}\,M_\odot$ and SFR $\approx10$–200 $M_\odot\,\mathrm{yr}^{-1}$, as derived from spatially integrated SED fitting with CIGALE, combining ground-based and Spitzer photometry from COSMOS2020 with JWST/NIRCam imaging from COSMOS-Web. We refer to Mitsuhashi et al. [69], Herrera-Camus et al. [43], and Faisst et al. [24] for full details on the SED-fitting procedure. High-resolution [C II] kinematic analyses from Lee et al. [58] are available for 16 out of 18 targets. The exceptions are DC–742174 (C–17) and VC–5110377875. A summary of the global properties of the galaxies is given in Appendix A.

JWST/NIRSpec IFU

The NIRSpec/IFU observation strategy and exposure times are outlined in Faisst et al. [24], while the data reduction pipeline used to derive science-quality data products is presented in Fujimoto et al. [34]. The survey adopts the medium resolution disperser/filter combinations G235M/F170LP ($\sim1.7$–$3.2\,\mu$m) and G395M/F290LP ($\sim2.9$–$5.3\,\mu$m), covering the main rest-frame optical emission lines (Hα, Hβ, [N II], [O II], [O III], and [S II]) at a spectral resolution of $R \approx1000$ ($\sim127$ km s$^{-1}$ in velocity dispersion). The data reduction pipeline largely follows the standard JWST/STScI pipeline (v1.16.0), as described in Rigby et al. [89], plus the application of several post-processing steps to improve the final data quality. Those routines consider additional cleaning of residual detector artifacts, improved background subtraction, masking of intermittent light leakages, removal of residual striping patterns in collapsed continuum maps, estimation of the data uncertainties, and astrometry alignment of the data cubes to the COSMOS-Web NIRCam mosaics, which are tied to the Gaia reference frame. We refer to Faisst et al. [24] and Fujimoto et al. [34] for more details.

DC–519281 / C–09

In the specific case of DC–519281 (or C–09), we analyze NIRSpec/IFU observations obtained under the Origin of the [C ii] Halos in Distant Systems (ORCHIDS) program (GO–5974; P.I.: M. Aravena). The main advantage of using ORCHIDS data is the improved spectral resolution due to the use of the G395H/F290LP disperser/filter combination ($R \approx2700$; $\sigma_{\mathrm{LSF}} \sim50\ \mathrm{km\ s}^{-1}$), plus the two-times finer pixel scale produced by the dithering strategy applied during the observations. The ORCHIDS data reduction follows the same pipeline used to reduce the ALPINE-CRISTAL-JWST sample data (Fujimoto et al. [34]), as briefly outlined above. We provide a comparison for both datasets in Appendix B.

JWST/NIRCam imaging

All targets are within the COSMOS field, benefiting from exquisite multi-wavelength coverage. High-quality JWST/NIRCam images (F115W, F150W, F277W, and F444W filters) from COSMOS-Web (PID 1727, PI: J. Kartaltepe, Casey et al. [9]) are available for most of our targets (17/18). A subset of the sample is also covered by PRIMER (PID 1837, PI: J. Dunlop, Donnan et al. [18]), providing additional imaging coverage (F090W, F200W, F356W, and F410M filters). At $z = 4 - 6$, these JWST images probe the rest-frame UV/optical continuum. We use F444W image mosaics drizzled onto a common grid with a pixel scale of 0′′.03. The F444W images trace the underlying stellar continuum (rest-frame $I$−band) of our targets. Full details on the JWST/NIRCam data reduction for the ALPINE-CRISTAL-JWST sample are given in Faisst et al. [24].

Methods

Our main goal is to study the Tully-Fisher relation, so we focus on carefully estimating the circular velocities ($V_{\mathrm{circ}}$) while adopting the stellar masses provided by Faisst et al. [24] for our targets (see § 2). We fit the data cubes pixel-by-pixel to build maps of H$\alpha$ emission-line intensity, line-of-sight velocity ($V_{\mathrm{LOS}}$), and velocity dispersion ($\sigma_V$). The kinematic data are then complemented with JWST F444W images to identify disk-like systems. For disk galaxies, we compute asymmetric drift correction terms and estimate $V_{\mathrm{circ}}$. Finally, to assess the effect of observational systematics we adopt a forward modeling approach.

Imaging analysis

We model the JWST/NIRCam F444W images to characterize the stellar morphology of our targets, as those roughly correspond to rest-frame $I$−band data for galaxies at $z = 4 - 6$ (Roper et al. [96]). We adopt a single-component Sérsic profile (Sérsic [98]) and use pysersic (Pasha & Miller [76]). The image model has seven free parameters: $x/y$ centroids, photometric position angle ($\mathrm{PA}_{\mathrm{phot}}$), effective radius ($R_e$), axis ratio ($q = b/a$), total flux, and Sérsic index ($n$). The fitting routine considers the convolution by the JWST/NIRCam point-spread function (PSF) taken from the JWST’s PSF simulation tool.¹ We note that this choice does not significantly affect the inferred parameters, as previous JWST/NIRCam studies have shown that simulated and

empirical PSFs lead to consistent results (Lee et al. [57]). Prior to modeling, we mask spurious features (e.g., foreground stars and corner artifacts) manually. When fitting, we found that the Sérsic indices were not well constrained due to the compactness of the sources and/or low S/N, but the best-fit values were typically close to unity (average $n = 1.03$, median of 0.95, and scatter of 0.47). Therefore, we fixed $n = 1$ (exponential profile). Fixing $n = 1$ produces only small changes in the morphological parameters, with median absolute differences of 7.2%, 1.6%, and 0.9% in $R_e$, inclination, and $\mathrm{PA}_{\mathrm{phot}}$, respectively. We show the quality of our Sérsic models in Appendix C.²

We assume a galactic oblate geometry to compute the inclination angles ($i$):

where $q$ is the projected axis ratio ($b/a$) of the galaxy on the sky, and $q_0$ is the intrinsic axis ratio of the galaxy. We assume $q_0 = 0.2$, a value appropriate for describing moderately thick disks (see e.g., Danhaive et al. [12]). This assumption is further motivated by observational evidence that high-$z$ star-forming galaxies have thick disks (van der Wel et al. [114]), as well as by cosmological simulations (e.g., Pillepich et al. [79]). Nevertheless, inferred kinematic quantities depend weakly on the exact choice of $q_0$ within commonly adopted values. For example, if the true intrinsic thickness were higher ($q_0 \approx0.25 - 0.5$), the inclination-corrected velocities would only vary by 10% (Danhaive et al. [12]), keeping the main results of this work unchanged (§ 4).

Emission line fitting

We adopt the H$\alpha$ emission line as our main tracer of the ionized-gas kinematics, as it enables a fair comparison between our measurements and pre-JWST TFR studies. For each datacube, we fit pixel-by-pixel the H$\alpha$, [N ii]$\lambda$6548, 6584, [S ii]$\lambda$6717, 6731 emission lines with Gaussian emission models using the LMFIT package (Newville et al. [73]) in Python. The underlying stellar continuum is modeled with a linear function as the observations’ S/N is found to be insufficient to detect any stellar absorption feature. The emission lines are fitted with single Gaussian components, tying their velocity centroid and velocity dispersion while leaving the individual line intensities free. We assume all lines trace the same bulk ionized-gas kinematics (e.g., Ho et al. (2016)). Hence, possible line-dependent centroid offsets or differences in linewidth are not captured by our model. Observed line widths are corrected for instrumental line broadening by considering the NIRSpec line-spread function (LSF) in convolution within the model. We adopt the LSF width measurements of Jakobsen et al. [46] for the grating/filter G395M/F290LP ($\sim127\ \mathrm{km\ s}^{-1}$ on average). To select only the spectra that are well-characterized by the emission line models, we compare the best-fit solutions with a simple straight-line (featureless) model and compute the Bayesian Information Criterion difference ($\Delta\mathrm{BIC}$, Liddle [62]). We keep pixels with $\Delta\mathrm{BIC} \equiv\mathrm{BIC}_{\mathrm{featureless}}-\mathrm{BIC}_{\mathrm{Gaussian}} > 10$ corresponding to “decisive evidence” in favor of the Gaussian emission-line model over the featureless model (Liddle [62]). After this step, we remove isolated spurious pixel detections by applying the binary_erosion task of scipy.ndimage³ with a $3 \times3$ pixel-box structure. However, as this task also erases good quality pixels characterizing high S/N emission from the galaxy boundaries, we recover those by applying the binary_dilation task of scipy.ndimage with the same pixel-box structure. The best-fit parameters are then used to build the Hα emission-line intensity, $V_{\rm LOS}$, and $\sigma_V$ two-dimensional maps (Figure 1). Only spaxels satisfying $\Delta\mathrm{BIC} > 10$ and retained by the morphological mask described above are included in the two-dimensional maps. No additional S/N threshold is applied.

Vector (PDF)

Disk classification

Given the relatively small sample size and the complex morphologies expected for galaxies at $z = 4$–$6$, we adopt a two-step classification scheme to identify disk-like galaxies in our sample. First, we visually inspect the Hα velocity fields and the JWST/NIRCam F444W images for all 18 galaxies, focusing on identifying close companions within the NIRSpec/IFU field-of-view (FoV; $5'' \approx30$–$35$ kpc at $z = 4$–$6$). We classify our targets into “galaxy pairs” and “isolated galaxy,” depending on whether a companion is present (or not) within the IFU field-of-view. We find that all the galaxy pairs show systemic velocity ($V_{\rm sys}$) difference of $\lesssim150$ km s$^{-1}$. Second, we apply the morpho-kinematic disk classification scheme of Wisnioski et al. [117] to the isolated galaxies. Briefly, the disk classification scheme of Wisnioski et al. [117] corresponds to five criteria that a source must satisfy for being labeled as “disk galaxy”, namely:

1. The velocity map exhibits a continuous velocity gradient along a single axis.

2. The ratio between the rotational velocity (amplitude of the velocity gradient) and the velocity dispersion must be greater than or similar to one.

3. The velocity dispersion map must peak near ($\lesssim1$ PSF FWHM) the kinematic center.

4. For inclined galaxies ($q < 0.6$), the photometric and kinematic position angles must agree within $< 30^\circ$.

5. The position of the kinematic center coincides, within one PSF, with the centroid of the continuum.

Photometric quantities are measured from F444W images (§ 3.1), while kinematic position angles (PA$_{\rm kin}$) are derived from H$\alpha$ V$_{\rm LOS}$ maps (§ 3.4). We support the assessment of criterion 3 by considering available [C ii] velocity dispersion profiles from Lee et al. [58]. The [C ii] line traces the cold, neutral gas component (see e.g. Algera et al. [1]); hence it provides an independent measurement of the velocity dispersion profile for the dominant mass component within our systems. We show the [C ii] velocity dispersion profiles as orange points in Figure 1. These profiles are extracted following the DYSMALpy-based procedure along PA$_{\rm kin}$ (see Lee et al. [58] for more details).

As a first step, we visually inspect the NIRSpec/IFU data cubes and the NIRCam/F444W images. We identify 14 isolated galaxies within the FoV of the IFU ($5'' \approx30 - 35$ kpc at $z = 4-6$). The other four systems correspond to galaxy pairs with $\Delta V_{\rm sys} \lesssim150$ km s$^{-1}$. Among the isolated galaxies, nine satisfy our first two disk-classification criteria, namely, they show a continuous velocity gradient aligned with a single axis and present a velocity amplitude-to-velocity dispersion ratio $\gtrsim1$. The third criterion requires the velocity dispersion to peak near the region of steepest velocity gradient, which is expected to be close to the kinematic center for disk-like systems. Only for DC–630594 we observe that the velocity dispersion peaks at the kinematic center, while the rest of the galaxies show velocity dispersion radial profiles that can be roughly regarded as constant. In practice, the ionized-gas velocity dispersion can be affected by non-circular motions (e.g., outflows) and by local features (Förster Schreiber et al. [32]; Phillips et al. [78]), making the interpretation of the maps complex. Thus, to provide an additional assessment, we complement the analysis of the H$\alpha$ kinematics with an examination of the [C ii] velocity dispersion profiles presented in Lee et al. [58]. The [C ii] line traces the total cold gas (Dessauges-Zavadsky et al. [14]) and, although it could also be affected by outflows (Ginolfi et al. [38]; Romano et al. [95]; Birkin et al. [6]), combining ionized gas and cold gas ($T \lesssim100$ K) information provides a more robust assessment of the dispersion peak. Using this approach, we find that seven galaxies satisfy the third criterion. The fourth criterion requires, for galaxies with $q < 0.6$, consistency between the photometric and kinematic major axes (within $30^\circ$), which is fulfilled by six galaxies. Finally, we discard the system VC–5100994794 because it does not meet the fifth criterion. That is, its kinematic center does not coincide with the continuum peak seen in the F444W image. Therefore, we find that a total of five galaxies can be regarded as disks ($5/18 = 27.8\%$), nine systems show perturbed kinematics ($9/18 = 50.0\%$), and the remaining four targets ($4/18 = 22.2\%$) correspond to galaxy pairs. Figure 2 shows our sample in the SFR–$M_\star$ plane separated by their morpho-kinematic classification (Table 1).

ID	Criteria	1	2	3	4$^{a}$	5	This work	[C ii]$^{b}$
DC–494763		✓	✓	✓	✓	✓	Disk	Best-Disk
DC–519281		✓	✓	✓	✓	✓	Disk	Disk
DC–536534		✓	✓	✓	✓	✓	Disk	Best Disk
DC–630594		✓	✓	✓	✓	✓	Disk	Best Disk
DC–742174$^{c}$		✓	✓	✓	✓	✓	Disk	…
DC–709575		✓	✓	×	…	…	Non-Disk	Non-Disk
DC–842313		×	…	…	…	…	Non-Disk	Non-Disk
DC–848185		×	…	…	…	…	Non-Disk	Disk
DC–873756		×	…	…	…	…	Non-Disk	Non-Disk
VC–5100541407		✓	✓	×	…	…	Non-Disk	Non-Disk
VC–5100994794		✓	✓	✓	✓	×	Non-Disk	Non-Disk
VC–5101218326		×	…	…	…	…	Non-Disk	Non-Disk
VC–5101244936		×	…	…	…	…	Non-Disk	…
VC–5110377875		✓	✓	✓	×	…	Non-Disk	…
DC–417567		×	…	…	…	…	Pair	…
DC–683613		×	…	…	…	…	Pair	Disk
DC–873321		×	…	…	…	…	Pair	Disk
VC–5100822602		×	…	…	…	…	Pair	Non-Disk

Vector (PDF)

We note that two of the galaxies classified as disk-like in this work have been previously discussed as possible AGN and/or ionized-gas outflow hosts. In particular, DC–536534 was identified as a robust AGN candidate by Ren et al. [87]. This target shows a broad H$\alpha$ emission line, regarded as a broad-line region component, plus a possible, additional outflow component. Ren et al. [87] also classify DC–519281 as an AGN candidate, with possible evidence for an outflow component bright in [O iii] emission. However, neither of these disk galaxies show significant cold-gas outflow signatures (Birkin et al. [6]).

We also note that DC–742174 is a borderline case in our disk classification scheme. On one hand, its rest-frame optical morphology shows two clump-like structures with no clear central emission. On the other hand, the H$\alpha$ velocity field displays a coherent velocity gradient along a single axis, the H$\alpha$ emission peaks in the galaxy center, and the observed velocity dispersion profile is broadly consistent with a roughly constant intrinsic dispersion. It is well-known that, at high-redshift, systems with irregular, often clumpy, morphologies in the rest-frame optical, can display regular velocity fields, reflecting an underlying disk nature (e.g., Förster Schreiber & Wuyts 2020, and references therein). However, given the limited spatial resolution of our observations, we cannot exclude the possibility that the observed kinematics arise from an unresolved interaction or from two systems projected close to each other in the sky. Therefore we flag DC–742174 separately throughout the analysis. Nevertheless, excluding DC–742174 from our analysis does not significantly change our results and conclusions (§ 4).

Kinematic measurements

For each disk galaxy, we start by characterizing its kinematics by deriving the major kinematic axis of the V$_{\rm LOS}$ map. The major kinematic axis is set by PA$_{\rm kin}$, which we determine by using PAFIT (Krajnović et al. [53]). PAFIT determines PA$_{\rm kin}$ by finding the axis that best reproduces the observed velocity field symmetry (see Krajnović et al. [53] for more details). Here, we fix the kinematic center to the centroid of the H$\alpha$ emission, as these are broadly consistent with the F444W photometric centroid within a PSF-sized region (full-width at half-maximum \sim0\farcs15-0\farcs17, D’Eugenio et al. [15]).

To construct rotation curves and $\sigma_V$ radial profiles, we follow the method used in DYSMALpy (Price et al. [82]; Lee et al. [59]). Briefly, we use the kinematic center and $\mathrm{PA}_{\mathrm{kin}}$ to derive the major kinematic axis. Across the major kinematic axis, we extract a set of integrated spectra using a series of circular apertures with a diameter of 0\farcs2 and separated by 0\farcs2, avoiding overlapping apertures and reducing covariance between neighbouring measurements. For each integrated spectrum, we fit the H$\alpha$ emission line to obtain its centroid and width, corresponding to the $V_{\mathrm{LOS}}$ and $\sigma_V$ values that we use to construct the rotation curve and the velocity dispersion radial profile (Figure 1). Here, we only consider best-fit results if the H$\alpha$ emission line has S/N $>5$.

From each rotation curve, we compute the inclination-corrected rotation velocity values by taking the average between $V_{\mathrm{LOS}}$ extracted at a radius $R=2.2R_{\mathrm{d}}$, where $R_{\mathrm{d}}$ is the disk-scale length.⁴ This choice is motivated by the fact that a self-gravitating exponential disk reaches its maximum rotation velocity at $\sim2.2R_{\mathrm{d}}$ (Freeman [33]). We note that this rotation velocity definition (often denoted $V_{2.2}$) is commonly adopted in H$\alpha$-based Tully-Fisher studies (e.g., Courteau [11]; Reyes et al. [88]).

We measure the global velocity dispersion of each galaxy from the corresponding aperture spectrum, correcting the observed line width for instrumental broadening using the LSF. We adopt an extraction radius of $2.2R_{\mathrm{d}}$ as a compromise between S/N and beam-smearing effects (e.g., Wisnioski et al. [117]; Johnson et al. [47]; Puglisi et al. [85]). That is, we avoid extracting values from the central regions where velocity gradients are strongest. We present the kinematic parameters for our galaxy disks in Table 2.

ID	$z_{\mathrm{H}\alpha}$	$R_{\mathrm{e}}$ (″)	$R_{\mathrm{e}}$ (kpc)	$i$ (°)	$\mathrm{PA}_{\mathrm{kin}}$ (°)	$V_{2.2}$ (km s$^{-1}$)	$V_{2.2}^{\mathrm{DYSMALpy}}$ (km s$^{-1}$)	$\sigma_0$ (km s$^{-1}$)	$\sigma_0^{\mathrm{DYSMALpy}}$ (km s$^{-1}$)	$V_{\mathrm{circ}}$ (km s$^{-1}$)	$V_{\mathrm{circ}}^{\mathrm{DYSMALpy}}$ (km s$^{-1}$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
DC–494763	5.2337	0.19	1.17	$53 \pm2$	$108 \pm7$	$109 \pm4$	$127 \pm4$	$51 \pm20$	$66 \pm13$	$153 \pm29$	$186 \pm20$
DC–519281	5.5759	0.16	0.95	$45 \pm3$	$319 \pm9$	$97 \pm5$	$113 \pm5$	$126 \pm39$	$129 \pm26$	$282 \pm77$	$315 \pm55$
DC–536534	5.6886	0.21	1.23	$57 \pm2$	$333 \pm3$	$80 \pm4$	$113 \pm3$	$115 \pm39$	$113 \pm23$	$254 \pm78$	$219 \pm32$
DC–630594	4.4403	0.22	1.46	$65 \pm2$	$143 \pm6$	$65 \pm6$	$110 \pm9$	$80 \pm29$	$78 \pm16$	$179 \pm56$	$195 \pm27$
DC–742174	5.6360	0.40	2.36	$76 \pm7$	$87 \pm10$	$90 \pm27$	$91 \pm3$	$84 \pm46$	$82 \pm16$	$197 \pm105$	$197 \pm31$

Asymmetric drift correction and circular velocity

Rotational velocities measured from ionized gas kinematics can be affected by pressure support, especially in high-$z$ galaxies where intrinsic velocity dispersions are high (Förster Schreiber & Wuyts 2020). This effect reduces the gas rotation velocity, so an asymmetric drift correction must be applied to $V_{2.2}$ to estimate $V_{\mathrm{circ}}$.

We follow the approach of Burkert et al. [8], where the circular velocity can be estimated as

where $R$ is the radius at which $V_{\mathrm{rot}}$ is derived. This asymmetric drift correction assumes: (i) a constant intrinsic velocity dispersion across the disk, (ii) an exponential surface mass density profile, and (iii) isotropic turbulent pressure support. We list the $V_{\mathrm{circ}}$ values in Table 2.

Forward modeling with DYSMALpy

The limited spectral and spatial resolution of the kinematic data of high-$z$ galaxies imposes severe limitations when recovering their intrinsic kinematic properties. Different kinematic routines have been developed to overcome such limitations, but systematic uncertainties inherent in the assumptions made by these tools can still bias the kinematic measurements (e.g., Lee et al. [58]; Yttergren et al. [119]). To address the effects of observational systematics on our kinematic parameters, we use the forward-modeling code DYSMALpy (Price et al. [82]; Lee et al. [59]) to derive the intrinsic kinematic parameters. The main advantage of DYSMALpy is that projection effects and observational broadening from the PSF and LSF are included directly in the forward modeling.

DYSMALpy fits composite mass models by generating mock kinematic observables, which can be a data cube, two-dimensional maps, or a rotation curve and a velocity dispersion profile, and compares them to the observed data. In this work, we fit the one-dimensional rotation curve and velocity dispersion radial profiles derived in § 3.4. The composite mass model can include a bulge, a disk, and a dark-matter halo. The baryonic mass component, bulge and disk, corresponds to the sum of two Sérsic profiles, where the bulge prominence is parameterized by the bulge-to-total mass ratio ($B/T$). However, since prominent bulges are not usual in typical star-forming galaxies at $z=4-6$ (Ferriera et al. 2023; Huerats-Company et al. 2024), and we do not observe any clear stellar features that could be regarded as bulges in the JWST/NIRCam F444W images, we set $B/T=0$. We also set the intrinsic thickness of the disk equal to $0.2$. We further assume a radially uniform intrinsic velocity dispersion, parameterized by a single value $\sigma_0$. For the dark-matter halo, we assume a two-parameter Navarro-Frenk-White profile (Navarro et al. [70]). We parameterize the halo contribution through the dark matter fraction at the effective radius, defined as $f_{\mathrm{DM}}(R_{\mathrm{e}})=V_{\mathrm{DM}}^2(R_{\mathrm{e}})/V_{\mathrm{circ}}^2(R_{\mathrm{e}})$. For each trial value of $f_{\mathrm{DM}}(R_{\mathrm{e}})$, DYSMALpy determines the corresponding NFW virial mass given the baryonic mass distribution and the adopted halo concentration. For each galaxy, we use the stellar-to-halo mass relation of Girelli et al. [39] to define a broad halo-mass range consistent with its stellar mass, adopting $M_{\mathrm{vir}}=10^{11}-10^{13}\,M_\odot$. We fix the concentration parameter ($c$) to 3.5, following the $M_{\mathrm{vir}}-c$ relation of Dutton & Maccio [20]. We compute the circular velocity by adopting the asymmetric drift approximation of Noordermeer [74] within DYSMALpy, so that the correction factor is calculated in the same way as described in Equation (2). The results obtained with DYSMALpy are presented in Table 2.

Results

Kinematic properties of the ALPINE-CRISTAL-JWST sample

Figure 1 presents the H$\alpha$ moment maps, rotation curves, and velocity dispersion radial profiles for the disk sample. The H$\alpha$ moment maps for the full sample are shown in Appendix D. The $V_{\mathrm{LOS}}$ maps show the typical blue-red pattern that indicates rotation. In all cases, the rotation curves and velocity dispersion radial profiles extend beyond $2.2R_{\mathrm{d}}$, allowing us to estimate $V_{2.2}$ and $\sigma_0$. For four galaxies, our inferred $\sigma_0$ estimates are lower than the instrumental resolution value ($\sim127$ km s$^{-1}$). This is because our emission-line modeling includes the LSF (§ 3) and the fitted emission lines have S/N high enough so that they are detected over several consecutive channels. We verify this by fitting mock spectra specifically designed to mimic the NIRSpec/IFU medium-filter spectral resolution. We find that intrinsic $\sigma_0$ values $\sim60-80$ km s$^{-1}$ can be recovered with typical errors below 40% if the fitted line has a S/N $\gtrsim10$ (see Appendix E

For more details. For the non-parametric measurements, the uncertainties inferred from these tests are included in the $\sigma_0$ error budget and propagated to the corresponding circular velocities.

The non-parametric measured kinematics described in § 3.4 can still be affected by beam smearing, even at the relatively large radii considered here. We therefore compare the non-parametric measurements with the intrinsic kinematic parameters inferred from the DYSMALpy forward models, which explicitly account for projection effects and convolution with the spatial PSF and spectral LSF. The intrinsic $V_{2.2}$ values inferred with DYSMALpy are systematically larger than the directly measured inclination-corrected values, by a factor of $\sim1.3$ on average ($\sim0.10$ dex). This shows that the limited spatial resolution suppresses the observed rotation velocities. In contrast, the intrinsic $\sigma_0$ values differ from the non-parametric measurements by only $\sim0.02$ dex on average, consistent within their uncertainties.

As independent checks, we also applied the empirical beam-smearing corrections of Johnson et al. [47] and performed independent 3D forward modelling with 3DBarolo (Di Teodoro & Fraternali [16]). Both methods imply smaller corrections to $V_{2.2}$, typically by a factor of $\sim1.1$, compared with the $\sim1.3$ correction inferred with DYSMALpy. These tests show that the observed rotation velocities are affected by the limited spatial resolution.

For the remainder of the analysis, we adopt the intrinsic kinematic parameters inferred with DYSMALpy as our fiducial values. Importantly, the differences among these beam-smearing treatments do not significantly affect our main results as discussed below.

Using our fiducial DYSMALpy measurements, the disk-like sample spans $V_{2.2} = 91 - 127$ km s$^{-1}$ and $\sigma_0 = 66 - 129$ km s$^{-1}$, corresponding to $V_{2.2}/\sigma_0 \simeq0.9 - 1.9$. These values indicate large random motions and significant asymmetric-drift corrections to be applied for estimating $V_{\mathrm{circ}}$. We discuss the impact of the adopted asymmetric-drift prescription on our results in Appendix F. After accounting for pressure support, we infer $V_{\mathrm{circ} } = 186 - 315$ km s$^{-1}$ for the disk-like sample (Table 2).

The stellar-mass Tully-Fisher relation

Figure 3 presents our targets in the context of the stellar-mass Tully-Fisher relation (sTFR). As a reference, we also show the local sTFR of Reyes et al. [88] and the sTFR at $z \approx2.3$ taken from Übler et al. [112]. We find that our disk galaxies are clearly offset from the local relation and cosmic-noon sTFRs.

Vector (PDF)

To quantify the average offset of our measurements with respect to the low-$z$ sTFRs, we fit a linear function of the form

where $\alpha$ is the slope and $\beta$ is the zero-point. Due to the small sample size, we cannot robustly constrain the slope. Therefore, we set $\alpha= 3.6$, the local universe value Reyes et al. [88], and fit only for the zero-point $\beta$, as commonly done in high-$z$ TFR studies (e.g., Puech et al. [84]; Tiley et al. [108]; Übler et al. [112]; Sharma et al. [100]). We perform the fit using lmfit, obtaining $\beta= 1.43 \pm0.18$. This value corresponds to a zero-point offset of $-0.96 \pm0.18$ dex with respect to the local sTFR of Reyes et al. [88] and $-0.50 \pm0.18$ dex lower relative to the $z \approx2.3$ sTFR of Übler et al. [112]. Varying $q_0$ from 0.2 to 0.5 changes the sTFR zero-point offset by 0.03 dex, below its statistical uncertainty. A similar variation of the sTFR zero-point is found when removing DC–742174 from the analysis.

We searched the literature for systems at $z = 4 - 6$ with JWST/NIRSpec IFU data delivering H$\alpha$ kinematics. We find four rotating or disk-like galaxies: a low-mass clumpy galaxy at $z \approx6$ with $M_\star= 4.5^{+2.7}_{-1.1} \times10^8\,M_\odot$ and SFR $= 2.6^{+1.7}_{-1.1}\,M_\odot\,\mathrm{yr}^{-1}$ (Fujimoto et al. [35]); GN20, an AGN host galaxy at $z \approx4$ with $M_\star= 1.1 \times10^{11}\,M_\odot$ (Übler et al. [111]); a rotating galaxy at $z = 4.26$ with $M_\star\approx4 \times10^{10}\,M_\odot$ and SFR $= 59 \pm5\,M_\odot\,\mathrm{yr}^{-1}$ (Jones et al. [48]) and another AGN host galaxy at $z = 4.76$ with $M_\star= 10^{10.98}\,M_\odot$ (Parlanti et al. [75]). When we include these systems in the $M_\star- V_\mathrm{circ}$ plane and re-fit the relation, we find $\beta= 1.54 \pm0.09$, corresponding to a zero-point offset of $-0.85 \pm0.09$ dex with respect to the local sTFR of Reyes et al. [88] and $-0.39 \pm0.09$ dex relative to the $z \approx2.3$ sTFR of Übler et al. [112]. The fitted zero-point changes by only $\Delta\beta= +0.11$ dex relative to the fit to our sample alone, which is smaller than the $1\sigma$ uncertainty of the latter. By using PYMC, we estimate a sTFR intrinsic scatter of 0.29 dex. We note, this value is not particularly driven by the inclusion of any of our targets or the literature data.

The baryonic-mass Tully-Fisher relation

In the local Universe, gas-rich low-mass galaxies are outliers of the sTFR. When the gas component is included in the mass budget, these galaxies follow a relation in terms of baryonic mass ($M_\mathrm{bar} = M_\star+ M_\mathrm{gas}$), further extending the sTFR in what has been called the baryonic-mass TFR (bTFR; McGaugh et al. [66]; Lelli et al. [61]). Galaxies at high-$z$ are commonly gas rich (Tacconi et al. [106]), suggesting that studying the bTFR may be better to determine any offset of the TFR zero-point. To estimate the cold-gas mass for our galaxies, we take the [C II] luminosities presented for the ALPINE galaxies (Béthermin et al. [3]) and adopt the $L_\mathrm{[CII]}-M_\mathrm{gas}$ prescription of Zanella et al. [120], with a scatter of $\sim0.3$ dex. We include this scatter in quadrature in the $M_\mathrm{bar}$ uncertainty. On average, we find that our galaxies are gas rich, with a mean sample gas fraction ($f_\mathrm{gas} = M_\mathrm{gas}/M_\mathrm{bar}$) of $0.68 \pm0.10$. This average gas fraction is consistent with the expected evolution of gas content in main-sequence galaxies with $M_\star\approx10^{10}\,M_\odot$ at $z \approx5$ (Tacconi et al. [106]), and with the mean $f_\mathrm{gas}$ value of the ALPINE parent sample ($f_\mathrm{gas} = 0.63 \pm0.03$, Dessauges-Zavadsky et al. [14]). These [C II]-based gas masses are also consistent with the dust-based estimates of Lee et al. [58] and obtained using the metallicity-dependent [C II] prescription of Vallini et al. [113]⁷. Relative to the Zanella et al. [120] estimates, we find mean offsets of $-0.26$ dex and $-0.18$ dex for Lee et al. [58] and Vallini et al. [113], respectively. Both offsets are consistent within the adopted $\sim0.3$ dex uncertainty of the Zanella et al. [120] conversion.

Figure 3 shows the bTFR for our disk galaxy sample. As a reference, we also show the local universe bTFR of Lelli et al. [61] and the bTFR at $z \approx2.3$ taken from Übler et al. [112]. We find all our sources located below the local bTFR. Three sources closely follow the cosmic noon bTFR, while the other two systems are offset from this relation. By fitting a linear function, keeping the slope fixed to the local value ($\alpha= 3.75$; Lelli et al. [61]), we find a zero-point offset of $-0.57 \pm0.20$ dex relative to the local bTFR, and $-0.30 \pm0.20$ dex relative to the $z \approx2.3$ bTFR. In both cases, the zero-point offsets are smaller compared to those measured for the sTFR. Similarly to the sTFR case, the literature data broadly agree with our sample. When we include these data in the fit, we obtain $\beta= 1.50 \pm0.08$, corresponding to a zero-point offset of $-0.69 \pm0.08$ dex with respect to the local bTFR of Lelli et al. [61] and $-0.42 \pm0.08$ dex relative to the $z \approx2.3$ bTFR of Übler et al. [112]. The fitted zero-point changes by $\Delta\beta= -0.12$ dex relative to the fit to our sample alone, which is smaller than the $1\sigma$ uncertainty of the latter. We note that, varying $q_0$ from 0.2 to 0.5 changes the inferred bTFR zero-point offset by 0.03 dex, and removing DC–742174 from our analysis changes the inferred zero-point offset by 0.07 dex. In both cases, the bTFR zero-point variation is well below the $1\sigma$ uncertainty.

We also perform a fit to the bTFR considering cold-gas masses inferred from the metallicity-dependent [C ii] prescription of Vallini et al. [113] and from the dust-based molecular gas fractions of Lee et al. [58]. For a consistent comparison, we use the four galaxies for which gas masses are available with all three prescriptions and compute the zero-point variation relative to that of assuming the conversion factor of (Zanella et al. [120]). The zero-point varies by $\Delta\beta= -0.14 \pm0.22$ dex and $-0.12 \pm0.34$ dex when adopting the Vallini et al. [113] and Lee et al. [58] prescriptions, respectively. These variations imply that our inferred bTFR zero-point offset is not driven by the specific cold-gas mass prescription adopted. By using pymc, we estimate a bTFR intrinsic scatter of 0.28 dex.

Discussion

Evolution of the TFR with cosmic time

Our study provides the first glimpse into the properties of disk galaxies and the Tully-Fisher relation at $z \approx5$. Our findings suggest a negative variation of the TFR zero point with redshift for both baryonic and stellar TFRs, meaning that galaxies at $z \approx5$ appear to have higher circular velocities for a fixed stellar or baryonic mass (Figure 3). Comparisons with galaxy samples at lower redshifts are not straightforward as studies have reported both weak and significant evolution of the sTFR (e.g. Puech et al. [84]; Vergani et al. [115]; Tiley et al. [107]; Espejo Salcedo et al. [23]; Mancera Piña et al. [64]). For example, Übler et al. [112] suggested a negative zero-point evolution up to $z \approx0.9$ for both bTFR and sTFR, and no evolution of the sTFR from $z \approx0.9$ to 2.3, which contrasts with the negative evolution of the bTFR reported over the same redshift range. Conversely, Tiley et al. [107] suggest only a mild negative evolution of the sTFR at $z \approx1$, consistent with Miller et al. [67, 68] and the predictions of N-body simulations of baryonic disks growing in dark matter halos of Dutton et al. [21]. More recently, Sharma et al. [100] reported a subtle evolution of the zero points of the bTFR and sTFR up to $z \sim0.6$–2.5. These contrasting results likely reflect differences in sample selection and/or kinematic modeling strategies when estimating $V_{\rm circ}$.

Interpreting the cosmic evolution of the zero points of the bTFR and sTFR is not straightforward. In the framework of virialized dark matter halos, the cosmic evolution of the zero point of the bTFR reflects the effect of smaller DM halo sizes at high-$z$, if the baryonic-to-halo mass ratio ($m_d \equiv M_{\rm bar}/M_{\rm halo}$) is assumed to remain constant with redshift. Because DM halo masses scale inversely with the Hubble parameter, $H(z)$, a smooth evolution is expected. This contrasts with the variation of the zero point of the sTFR with redshift, which is more convoluted due to the cosmic evolution of the $f_{\rm gas}$ (Tacconi et al. [106]), which also plays a role (Übler et al. [112]). For both bTFR and sTFR, further considerations regarding the internal mass distribution of galaxies must be made because $V_{\rm circ}$ is usually estimated at inner radii (e.g., $\approx2.2R_d$), not at the virial radius of the halo. The low DM fraction measured in the central regions of galaxies at high-$z$ (Genzel et al. [36]; Lang et al. [55]; Genzel et al. [37]) biases the central $V_{\rm circ}$ to low values, shifting galaxies leftward in the $M_{\rm bar}$–$V_{\rm circ}$ plane, such that any possible evolution of the TFR could be erased or even reversed (Übler et al. [112]). To account for these effects, we compare our measurements with the predictions of the empirical model of Übler et al. [112]:

where $f_{\rm DM}(<R_{\rm e},z)$ is the dark matter fraction within one effective radius. The model assumes virialized DM halos, galaxy sizes proportional to their parent halo sizes, and baryonic masses related to those of the halos through $f_{\rm DM}$ (see Appendix D of Übler et al. [112] for its full derivation, and Mancera Piña et al. [64] for a similar model). In this toy model, the variation of the TFR zero point with redshift can be further affected by (i) changes in the baryonic-to-halo mass ratio, (ii) the evolution of the DM fraction in the inner regions of the halo (parameterized at the optical radius), and (iii) changes in the gas fraction for the sTFR. Initially, this model was introduced to study the evolution of the TFR at cosmic noon by Übler et al. [112]. Here, we update it by adopting more recent parameterizations of $f_{\rm DM}(R_{\rm e},z)$, $f_{\rm gas}(z)$, and $m_d(z)$ for main-sequence galaxies with $M_\star= 10^{10} M_\odot$ and $z \leq5$. We adopt the parameterization for $f_{\rm gas}(z)$ from Tacconi et al. [106]. To derive $m_d(z)$, we assume a stellar-to-halo mass ratio combined with $f_{\rm gas}(z)$. We assume a constant stellar-to-halo mass ratio of $0.020 \pm0.005$ due to its expected modest cosmic evolution for galaxies with $M_\star= 10^{10} M_\odot$ (Girelli et al. [39]). To compute $f_{\rm DM}(R_{\rm e},z)$, we combine measurements from the literature and interpolate them with redshift. We adopt $f_{\rm DM}(<R_{\rm e}) \approx0.64$ for galaxies at $z \lesssim0.6$ (Martinson et al. [65]). At $0.6 \lesssim z \lesssim2.5$, we consider $f_{\rm DM}(R_{\rm e},z)$ values taken from the RC100 sample (Nestor Shachar et al. [72]). Finally, for systems at $z \sim5$, we include the average estimate for CRISTAL galaxies (Lee et al. [58]). We interpolate these values by fitting a smooth redshift function, obtaining $f_{\rm DM}(z) = 0.58 \exp[-(0.31 \times z)^{1.5}]$ (see Appendix G for more details).

Figure 4 shows the zero-point variation with redshift for the bTFR and sTFR. As a reference, we also include the expected trends for virialized halos in the bTFR, scaling as the inverse of $H(z)$, and its conversion to the sTFR when accounting for the cosmic evolution of $f_{\rm gas}$. In both cases, these predicted trends overestimate the zero-point offsets measured by $\sim-0.3$ dex at $z \sim5$, unlike the empirical model, which remains consistent within $1\sigma$ uncertainty. The updated and the original empirical model of Übler et al. [112] predict different variations in the bTFR zero-point offset with redshift, highlighting how sensitive these are to the adopted $f_{\rm DM}(<R_{\rm e},z)$ prescription (Appendix G). However, in both model versions, the relatively small cosmic evolution of the bTFR (and sTFR) is due to high-$z$ massive galaxies presenting low DM fractions in their inner regions.

Vector (PDF)

Our empirical model does not include any variation in the intrinsic properties of the DM halos (e.g., concentration, spin) with redshift (Übler et al. [112]). Such variations could occur depending on how baryons are accreted and assembled in the inner halo regions (Glowacki et al. [40]; Sorini et al. [103]). For example, decreasing halo concentrations with redshift could reduce the cosmic evolution of the bTFR zero point, as highlighted by the semi-analytic model (SAM) of Dutton et al. [21]. Although that model predicts zero-point offsets only up to $z \approx4$ for the bTFR and sTFR, these trends broadly agree with the empirical model of Übler et al. [112] for $z \gtrsim2$ systems, and are also broadly consistent with our results.

We caution that the inferred zero-point offsets are based on a small number of disk galaxies. Although fixing the slope reduces the number of free parameters, small samples can still introduce statistical fluctuations in the recovered TFR normalization. To estimate this effect, we performed recovery tests by generating mock linear TFRs with a known input offset and re-fitting randomly selected subsamples of five and nine galaxies. On average, the input offsets were recovered without significant bias, with a mean difference between the recovered and input zero-points of $\sim0.01$ dex. The scatter in the recovered offsets was $\sim0.25$ dex for five galaxies and $\sim0.18$ dex for nine galaxies, both smaller than the offsets reported here. This suggests that the measured negative evolution is unlikely to be driven entirely by small-number statistics. Nevertheless, the exact amplitude of the evolution should be regarded as tentative until larger samples of spatially resolved disk galaxies become available.

Are we witnessing the settling of the first galactic disks?

Recent studies of galaxy morphology ([29, 50, 45]) and cold-gas kinematics ([58]) report that $\sim30\% - 50\%$ of the $z \approx5$ systems resemble disk-like structures. We also find a comparable disk fraction ($f_{\mathrm{disk}}$) in our sample, with five disk-like systems out of 18, corresponding to $f_{\mathrm{disk}} = 0.28 \pm0.11$.⁸ If we also consider galaxies in pairs which fulfill our disk criteria, we find a total of eight disk-like systems, yielding $f_{\mathrm{disk}} = 0.36 \pm0.10$ (Figure 5). These values indicate that a non-negligible subset of massive galaxies ($M_{\star} \gtrsim10^{9.5}$) has already developed a rotating component by $z \approx5$, while a larger fraction of the population still exhibits disturbed kinematics and/or is undergoing merging.

Vector (PDF)

It remains unclear how these disk galaxies arise in the early universe and the main mechanisms driving their evolution. On one hand, frequent interactions can perturb and disrupt forming disks, so observations of both isolated and interacting disks can be interpreted as galaxy growth occurring in environments where major encounters remain common ([92, 86, 19]). On the other hand, high gas accretion rates and strong internal gravitational instabilities can lead to increased turbulence within the disks, keeping them evolving in a gravitational state marginally stable (e.g., [117, 104]). The average velocity dispersion ($\sim90$ km s$^{-1}$) and $V_{2.2}/\sigma_{0}$ ratio ($= 1.26 \pm0.19$) measured for our disk galaxies are broadly consistent with the latter scenario (Figure 6), as they approximately follow the cosmic evolutionary trend predicted by the marginally stable disk model of Wisnioski et al. [117], which assumes a Toomre (1964) parameter $Q \approx1$ (see Bacchini et al. [2] for some limitations of this model). Even though we do not directly compute $Q$, the rough agreement with the cosmic evolutionary trend suggest that the kinematics of these galactic disks at $z \approx5$ likely reflect the balance between rotation, pressure support, and feedback-driven turbulence as they grow and settle while their gas supply changes ([51, 117, 102]).

Vector (PDF)

The kinematic tracer we used, Hα, could include some contribution from non-circular motions, such as ionized gas outflow or extraplanar regions. As a result, intrinsic velocity dispersion from ionized gas may overestimate the disk mid-plane value, which can be more reliably recovered from cold-gas kinematics (Wisnioski et al. [118]). However, comparisons with available [C II] kinematic data from the CRISTAL survey suggest that this is not the case. For the four disk-like galaxies in common with Lee et al. [58], we find that their average intrinsic cold-gas velocity dispersion is $\sim87$ km s$^{-1}$, consistent with our Hα-based measurements ($\sim93$ km s$^{-1}$). The full CRISTAL disk sample has a median intrinsic dispersion of $\sim70$ km s$^{-1}$, consistent with the high gas fractions of systems and turbulence being driven by gravitational instabilities and gas transport (Lee et al. [58]). This is somewhat different from the predictions of post-processed simulations (Kohandel et al. [52]), which suggest a factor of $\sim2$ difference between the velocity dispersion of the ionized and molecular gas ISM phases for the stellar-mass range probed by our disk galaxy sample ($M_\star\approx10^{9.5-10.5}\,M_\odot$). More generally, the empirical evidence for large tracer-dependent offsets in disk kinematics for high redshift galaxies remains inconclusive. The most direct test is to compare multiple tracers within the same galaxies. When this has been done, no clear systematic offset in $V_{\mathrm{rot}}/\sigma_0$ between different gas phases has been found ([58, 118]).

It remains unclear which ISM phase predominantly gives rise to the observed [C II] emission, and whether the relative contributions from different ISM phases vary with underlying galaxy properties (Wisnioski et al. [118]). Another possibility is that

the [C II]-based velocity dispersions are overestimated due to the modest spatial resolution of the interferometric data ($\sim2$–4 kpc; see Lee et al. [58] for more details).

Overall, the galaxy disks identified here display ordered rotation with moderate values of $V_{2.2}/\sigma_0 \approx1$, indicative of dynamically warm systems. We note that these estimates are consistent with the JWST/slitless (grism) measurements made for a larger population of $\sim100$ galaxies at $z = 4$–6 (Danhaive et al. [12], $\sigma_0 \approx100$ km s$^{-1}$ and $V_{\mathrm{rot}}/\sigma_0 \approx1.0$). Our measurements suggest that these early disks are not yet thin, cold disks like those we commonly observe at low redshift (Epinat et al. [22]). Instead, they appear thick and turbulent, likely sustained by rapid gas accretion and intense feedback events (Krumholz et al. [54]). The fact that we detect a significant disk fraction at $z = 4$–6, and that their kinematics can be understood under the framework of gravitationally marginally stable turbulent disks, provides clues that disk settling is already underway. Previous works have already suggested that massive galaxies ($M_\star\gtrsim10^{10}\,M_\odot$) at $z \approx4$–5 can already be settled on rotating disks ([71, 60]). In contrast, the disk candidates identified in this work are less massive, with $M_\star\gtrsim10^{9.5}\,M_\odot$. The broad diversity of kinematics and morphologies seen in the full ALPINE-CRISTAL-JWST sample suggests that disk settling is an emerging process at this cosmic epoch and over a wider range of stellar mass than previously thought.

Conclusions

We have presented a first analysis of the stellar mass and baryonic Tully-Fisher relations at $z = 4$–6 using JWST/NIRSpec IFU observations from the ALPINE-CRISTAL-JWST survey, complemented with JWST/NIRCam imaging. Our main results can be summarized as follows:

1. We identify five out of 18 galaxies as isolated disk-like systems, corresponding to a sample disk fraction of $f_{\mathrm{disk}} = 0.28 \pm0.11$ at $z = 4$–6. In addition, four systems are classified as galaxy pairs, while the remaining nine show perturbed kinematics.

2. The five disk-like galaxies show clear velocity gradients and rotation curves that extend up to or beyond $\sim2.2R_{\mathrm{d}}$, allowing estimates of their kinematic properties. The sample spans intrinsic $V_{2.2}=91-127$ km s$^{-1}$, $\sigma_0=66-129$ km s$^{-1}$, and $V_{2.2}/\sigma_0=0.9-1.9$, suggesting that random motions remain important in these systems. Our measured intrinsic velocity dispersions and $V_{2.2}/\sigma_0$ ratios are consistent with the expectations for marginally stable disks given the average gas fraction of the systems at this cosmic epoch.

3. We find that the disk galaxies in our sample are offset with respect to the low-$z$ stellar-mass and baryonic Tully-Fisher relations. In the sTFR, all five disks lie below the local relation and the one at $z\approx2.3$, suggesting zero-point offsets of $-0.96\pm0.18$ dex and $-0.50\pm0.18$ dex, respectively. In the bTFR, all five disks lie below the local relation, while three are broadly consistent with the $z\approx2.3$ bTFR. We estimate zero-point offsets of $-0.57\pm0.20$ dex relative to the local bTFR and $-0.30\pm0.20$ dex relative to the $z\approx2.3$ relation.

4. We extend the empirical model of Übler et al. [112] up to $z=5$, finding good agreement with our best-fit zero-points of both the sTFR and bTFR. In this framework, the redshift evolution of both TFR zero-points is driven by the cosmic evolution of dark matter halo properties, modulated by a reduced contribution of dark matter to the dynamical mass on the galaxy scale.

Overall, our results suggest that rotating disk-like systems were already in place within the first $\sim1$ Gyr of cosmic time, although their kinematics show high degree of turbulence. While our conclusions are limited by small number statistics, this work provides one of the first spatially resolved views of the Tully-Fisher relation for main-sequence galaxies at $z=4-6$. Larger JWST/NIRSpec IFU samples, together with deeper ALMA observations for the cold gas and spatially resolved dynamical modeling, will be key for determining whether the offsets reported here are representative of the bulk galaxy population at these redshifts.

Acknowledgements. We thank the anonymous referee for their careful reading and constructive comments, which greatly improved the quality and clarity of this manuscript. We warmly thank Annagrazia Puglisi and Andrés Escala for their comments, which helped to improve this study. D.G.E. acknowledges support from ANID BECAS/DOCTORADO NACIONAL/2024-21251071. D.G.E., J.M., E.I., M.A., R.H.C and V.V. gratefully acknowledge support from ANID MILENIO NCN2024_112, D.G.E. and E.I. acknowledge support from ANID FONDECYT Regular 1221846. M.A. is supported by FONDECYT grant number 1252054. M.A. acknowledges support from ANID+Vinculación Internacional + FOVI250261. M.A., B.B., RJA, and V.V. acknowledge support from the ANID BASAL project FB210003. This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence of Université Côte d’Azur under reference number ANR-15-IDEX-01. This research was funded, in whole or in part, by the French National Research Agency (ANR), grant ANR-25-CE31-0192 (project CIGALE2D) RJA was supported by FONDECYT grant number 1231718. EdC acknowledges support from the Australian Research Council through project DP240100589. M. Palla acknowledges support from HORIZON-INFRA-2024-DEV-01-01 – Research Infrastructure Concept Development, through the project WST: The Wide-Field Spectroscopic Telescope (Grant No. 101183153). M. Relano acknowledges support from project PID2023-150178NB-I00 financed by MCIU/AEI/10.13039/501100011033, and by FEDER, UE. HÜ thanks the Max Planck Society for support through the Lise Meitner Excellence Program. HÜ acknowledges funding by the European Union (ERC APEX, 101164796). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

1. Instituto de Física y Astronomía, Universidad de Valparaíso, Avda. Gran Bretaña 1111, Valparaíso, Chile.

2. Millennium Nucleus for Galaxies (MINGAL).

3. Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago 8370191, Chile.

4. Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, 06000, Nice, France.

5. Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France.

6. Dipartimento di Fisica e Astronomia Galileo Galilei Università degli Studi di Padova, Vicolo dell’Osservatorio 3, 35122 Padova Italy.

7. INAF - Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122, Padova, Italy.

8. International Centre for Radio Astronomy Research (ICRAR), The University of Western Australia, M468, 35 Stirling Highway, Crawley, WA 6009, Australia.

9. Département d’Astronomie, Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland.

10. IPAC, California Institute of Technology, 1200 E. California Blvd. Pasadena, CA 91125, USA.

11. Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA.

12. David A. Dunlap Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada.

13. Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile, Santiago 7820436, Chile.

14. Observatory, Carnegie Institution of Washington, Raúl Bitrán 1200, La Serena, Chile.

15. Departamento de Astronomía, Universidad de Concepción, Barrio Universitario, Concepción, Chile.

16. Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan.

17. Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA.

18. Max-Planck-Institute für extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany.

19. Department of Physics and Astronomy and George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, 4242.

20. Department of Physics and Astronomy, University of California, Riverside, 900 University Ave, Riverside, CA 92521, USA.

21. INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy.

22. Dept. Física Teórica y del Cosmos, Campus de Fuentenueva, Edificio Mecenas, Universidad de Granada, E-18071, Granada, Spain.

23. Instituto Universitario Carlos I de Física Teórica y Computacional, Universidad de Granada, 18071, Granada, Spain.

24. Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany.

25. Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI, UTIAS, Tokyo Institutes for Advanced Study, University of Tokyo, Chiba, 277-8583, Japan.

26. Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan.

27. Center for Data-Driven Discovery, Kavli IPMU (WPI, UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan.

28. Center for Astrophysical Sciences, Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA.

29. INAF – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti 93/3, 40129 Bologna, Italy.

Vector (PDF)

Vector (PDF)

Global properties of the ALPINE-CRISTAL-JWST sample

In Table A.1, we present a compilation of different properties for the full ALPINE-CRISTAL-JWST sample.

Comparison of datasets for CRISTAL-09

The available ORCHIDS observations for DC–519281/C–09 allow us to assess whether the ALPINE-CRISTAL-JWST data are sufficient for recovering the intrinsic velocity dispersion values for this source. In Figure B.1 we present a comparison between these two datasets. The ORCHIDS observations use G395H/F290LP ($R \approx2700$) with a spatial sampling of $0.05^{\prime\prime}$ pixel$^{-1}$, while the ALPINE-CRISTAL-JWST observations use G395M ($R \approx1000$) with $0.10^{\prime\prime}$ pixel$^{-1}$ sampling. Overall, we find broad agreement. The moment maps show consistent features among both datasets, despite the enhanced data quality delivered by the observing setup of ORCHIDS. For a consistent comparison, the kinematic profiles are extracted using circular apertures of $0.2^{\prime\prime}$ in diameter and separated by $0.2^{\prime\prime}$ along the same kinematic major axis. The aperture positions are defined from the G395H data with S/N $> 5$ and the same positions are used for G395M. One of the four G395H apertures has S/N $< 5$ in G395M and is therefore not included in the G395M profile.

This agreement can be quantitatively assessed in Figure B.2, where the rotation velocities agree within $1\sigma$ at two of the three common aperture positions, and the velocity dispersions independently agree within $1\sigma$ at two of the three positions. The $V_{\mathrm{los}}$ extracted at $2.2R_{\mathrm{d}}$ is $55 \pm3$ km s$^{-1}$ for the G395H grating and $55 \pm5$ km s$^{-1}$ for the G395M grating. The velocity dispersion at $2.2R_{\mathrm{d}}$ is $128 \pm3$ km s$^{-1}$ for G395H and $113 \pm7$ km s$^{-1}$ for G395M. The uncertainties tend to be smaller in the high resolution data.

Exponential profiles for the full sample

Figure C.1 shows the F444W images, the adopted exponential ($n = 1$) PySérsic models, and the corresponding residuals for the full sample.

Non-disks and galaxy pairs

For completeness, Figure D.1 shows the JWST false-colour images and Hα moment maps for the non-disk systems and galaxy pairs.

Recovery test for $\sigma_V$

Some of the velocity dispersion values for our sources are lower than the LSF width of the NIRSpec-IFU G395M grating ($\sigma_{\mathrm{LSF}} \approx127$ km s$^{-1}$). This is because many observed emission lines are detected over several consecutive spectral channels due to their high S/N. Nevertheless, it is important to determine the minimum S/N at which we can reliably recover $\sigma_V$ given the limited spectral resolution ($R \approx1000$).

To address this, we performed a set of simulations adapted to the properties of the data for each galaxy. First, we generated synthetic Hα emission lines, assuming Gaussian profiles and adopting different intrinsic velocity dispersion values. Then, the simulated emission lines were shifted to the redshift of each galaxy, convolved with the instrumental LSF (Shajib et al. [99]), and resampled to match the spectral step of the JWST data cube. To reproduce realistic observational conditions, we added noise extracted directly from emission-line-free spaxels in the data cubes. We produced many mock emission lines with varying intrinsic velocity dispersions and S/N ratios, considering multiple noise extraction realizations. Finally, these mock spectra were analyzed using the same procedure applied to the observed emission lines (§ 3). Figure E.1 compares the intrinsic velocity dispersion with the recovered values. The G395M grating data quality is insufficient to recover emission lines with $\sigma_V \lesssim40$ km s$^{-1}$, regardless of the observed S/N. In the regime most relevant for our data (S/N $\approx5$–20 and intrinsic velocity dispersions $\sim70$–100 km s$^{-1}$), we find that the velocity dispersion can be recovered with an accuracy $\sim20$–40%. This result indicates that we can reliably recover $\sigma_V$ for our systems and is consistent with the comparison between different JWST NIRSpec-IFU observing setups for DC–519281.

Asymmetric-drift correction

Given the comparable values of $V_{2.2}$ and $\sigma_0$ in our sample, asymmetric-drift corrections are significant: the circular velocities are obtained by applying factors of $\sim2$–3 to the rotational velocities. However, the magnitude of the asymmetric-drift corrections depends on estimates for the pressure gradient within galactic disks, which are based on assumptions about the gas surface density, velocity dispersion, and disk vertical structure (see appendix in Bouché et al. [7]). These estimates are difficult to constrain at high redshift.

ID	CRISTAL ID	R.A. (J2000.0)	Decl. (J2000.0)	$z_{\mathrm{H}\alpha}$	$\log M_\star$ ($M_\odot$)	$\log\mathrm{SFR}_{\mathrm{H}\alpha}$ ($M_\odot\,\mathrm{yr}^{-1}$)	$\log L_{[\mathrm{CII}]}$ ($L_\odot$)	$\log M_{\mathrm{gas}}$ ($M_\odot$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
DC–417567	C–10	10:02:04.08	+01:55:44.0	5.6700	9.99 ± 0.31	1.86 ± 0.02	8.50 ± 0.08	9.98 ± 0.32
DC–494763	C–19	10:00:05.10	+02:03:12.2	5.2337	9.51 ± 0.36	1.76 ± 0.03	8.69 ± 0.05	10.17 ± 0.31
DC–519281	C–09	09:59:00.89	+02:05:27.6	5.5759	9.84 ± 0.39	2.12 ± 0.05	8.73 ± 0.09	10.21 ± 0.32
DC–536534	C–03	09:59:53.24	+02:07:05.4	5.6886	10.40 ± 0.29	2.61 ± 0.04	8.90 ± 0.09	10.39 ± 0.32
DC–630594	C–11	10:00:32.60	+02:15:28.4	4.4403	9.68 ± 0.32	1.83 ± 0.08	8.79 ± 0.04	10.27 ± 0.31
DC–683613	C–05	10:00:09.42	+02:20:13.8	5.5420	10.16 ± 0.35	1.60 ± 0.05	8.90 ± 0.04	10.39 ± 0.31
DC–709575	C–14	09:59:47.06	+02:22:32.8	4.4121	9.53 ± 0.38	1.36 ± 0.06	8.48 ± 0.09	9.96 ± 0.32
DC–742174	C–17	10:00:39.12	+02:25:32.2	5.6360	9.51 ± 0.40	1.18 ± 0.06	8.17 ± 0.10	9.64 ± 0.32
DC–848185	C–02	10:00:21.51	+02:35:10.9	5.2931	10.30 ± 0.28	2.43 ± 0.03	9.20 ± 0.03	10.69 ± 0.31
DC–873321	C–07	10:00:04.06	+02:37:35.8	5.1542	10.00 ± 0.33	2.05 ± 0.03	8.98 ± 0.05	10.47 ± 0.31
DC–873756	C–24	10:00:02.72	+02:37:40.1	4.5457	10.53 ± 0.08	0.30 ± 0.50	9.56 ± 0.02	11.06 ± 0.31
DC–842313	C–01	10:00:54.53	+02:34:35.0	4.5537	10.65 ± 0.50	2.50 ± 0.14	8.66 ± 0.08	10.14 ± 0.32
VC–5100541407	C–06	10:01:00.92	+01:48:33.7	4.5630	10.09 ± 0.30	0.70 ± 0.49	9.07 ± 0.04	10.56 ± 0.31
VC–5100822662	C–04	09:58:57.91	+02:04:51.4	4.5205	10.15 ± 0.29	1.83 ± 0.08	8.90 ± 0.04	10.39 ± 0.31
VC–5100994794	C–13	10:00:41.14	+02:17:14.1	4.5802	9.65 ± 0.33	1.89 ± 0.13	8.75 ± 0.04	10.23 ± 0.31
VC–5101218326	C–25	10:01:12.50	+02:18:52.6	4.5739	10.90 ± 0.32	1.04 ± 0.49	9.26 ± 0.02	10.75 ± 0.31
VC–5101244936	C–15	10:00:47.65	+02:18:02.1	4.5803	9.69 ± 0.33	1.88 ± 0.09	8.70 ± 0.09	10.18 ± 0.32
VC–5110377875$^{(\ast)}$	…	10:01:32.34	+02:24:30.3	4.5505	10.17 ± 0.17	1.96 ± 0.17	9.23 ± 0.03	10.72 ± 0.31

Here, we test how different asymmetric-drift prescriptions affect the inferred stellar Tully-Fisher offset. We start by rewriting Equation (2) as

where $\alpha$ is the factor that varies depending on the adopted prescription. A galaxy model with a constant scale height gives $\alpha= 1$ (Bouché et al. [7]). For thick-disk turbulent galaxies with mass distributions following flattened Sérsic profiles, we have $\alpha\approx1.3$ for a self-consistent asymmetric-drift correction (Price et al. [83]). Finally, the Burkert et al. [8] prescription adopted in our main analysis gives $\alpha= 2$, which assumes an exponential gas surface density profile and a constant, isotropic velocity dispersion.

Fig. F.1 shows the stellar Tully-Fisher offset relative to Reyes et al. [88] for these three prescriptions. Lower values of $\alpha$ produce less negative offsets, as expected. However, all prescriptions agree within the uncertainties. We adopt the Burkert et al. [8] correction because it is consistent with our use of a single characteristic dispersion, $\sigma_0$, and allows a direct comparison with previous high-redshift Tully-Fisher studies.

Cosmic evolution of $f_{\mathrm{DM}}(< R_{\mathrm{e}})$

The redshift evolution of $f_{\mathrm{DM}}(< R_{\mathrm{e}}, z)$ is the least constrained variable in the model of Übler et al. [112]. Here, we derive a functional form for $f_{\mathrm{DM}}(< R_{\mathrm{e}}, z)$ by fitting measurements from the literature. At $z = 0$, we consider $f_{\mathrm{DM}} = 0.64 \pm0.18$ from the DiskMass Survey (Martinsson et al. [65]). For the $z = 0.7 - 2.5$ range, we adopt the measurements for the RC100 sample published in Nestor Shachar et al. [72]. It is worth noting that the stellar-mass ranges of those comparison samples are slightly higher than the range probed in this work, $10^{10-11}\,M_\odot$ versus $10^{9.5-10.5}\,M_\odot$. Nevertheless, we adopt these samples because they currently represent the lowest-mass systems available in the literature with reliable $f_{\mathrm{DM}}$ determinations at $z = 1 - 3$. Additionally, we take the average $f_{\mathrm{DM}} = 0.08 \pm0.07$ value derived for our disk galaxies in Lee et al. [58]. During the fitting, we adopt the following parameterization:

where $A$, $B$, and $C$ are free parameters. We find

We show the best-fit results in Figure G.1. By comparing with the power-law parameterization of Nestor Shachar et al. [72] for the RC100 sample, we find that our best-fit predicts lower $f_{\mathrm{DM}}(< R_{\mathrm{e}}, z)$ values at $z \gtrsim2.5$. At low-$z$, our best-fit also predicts lower dark matter fractions, but the $f_{\mathrm{DM}}(< R_{\mathrm{e}}, z)$ trend does not peak at high values as in the parameterizations of Nestor

Vector (PDF)

Shachar et al. [72]. The parameterization of Übler et al. [112] considers higher $f_{\mathrm{DM}}(< R_{\mathrm{e}})$ values for galaxies at low-$z$, but it broadly agrees with ours at $z \gtrsim2.3$.

In practice, the $f_{\mathrm{DM}}(< R_{\mathrm{e}})$ data show no preference for any model due to the significant uncertainties. However, compared with the parameterization of Nestor Shachar et al. [72] our best-fit $f_{\mathrm{DM}}(< R_{\mathrm{e}},z)$ induces a substantially different prediction for the bTFR (and sTFR) zero-point evolution with redshift in the model of Übler et al. [112]. The parameterization of Nestor Shachar et al. [72] implies positive offsets for the bTFR at low-$z$, while our best-fit does not. This is because the $f_{\mathrm{DM}}(< R_{\mathrm{e}})$ parameterization of Nestor Shachar et al. [72] accounts for a significant reduction in the dark matter content within galaxies at $z \lesssim1$, which counterbalances and further dominates the effect induced by the $H(z)^{-1}$ term in the model of Übler et al. [112], see Equation 4a. This effect can also be seen when adopting the $f_{\mathrm{DM}}(< R_{\mathrm{e}})$ parameterization of Übler et al. [112]. The positive zero-point evolution of the bTFR (and sTFR) described by assuming the $f_{\mathrm{DM}}(< R_{\mathrm{e}})$ parameterization Nestor Shachar et al. [72] is not consistent with the observational data at $z \lesssim1$, which is why we prefer our analytical prescription for $f_{\mathrm{DM}}(< R_{\mathrm{e}},z)$.

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)