11institutetext: Department of Physics and Astronomy, University of the
Western Cape, Cape Town 7535, South Africa 22institutetext: SISSA
International School for Advanced Studies, Via Bonomea 265, I-34136 Trieste,
Italy 33institutetext: QGSKY, INFN-Sezione di Trieste, via Valerio 2, I-34127
Trieste, Italy 44institutetext: IFPU Institute for Fundamental Physics of the
Universe, Via Beirut, 2, 34151 Trieste, Italy 55institutetext:
Universit$\acute{e}$ de Strasbourg, CNRS UMR 7550, Observatoire astronomique
de Strasbourg, 67000 Strasbourg, France 66institutetext: Department of
Astrophysics, University of Vienna, Türkenschanzstraße 17, 1180 Vienna,
Austria 77institutetext: Sterrenkundig Observatorium, Universiteit Gent,
Krijgslaan 281 S9, 9000 Gent, Belgium 88institutetext: Department of Astronomy
and Theoretical Physics, Lund Observatory, Box 43, SE-221 00 Lund, Sweden
99institutetext: University of Strasbourg Institute for Advanced Study, 5
allée du Général Rouvillois, F-67083 Strasbourg, France

#  Dark Matter Fraction in Disk-Like Galaxies Over the Past 10 Gyr

Gauri Sharma  , Contact<EMAIL_ADDRESS><EMAIL_ADDRESS>Jonathan
Freundlich 55 Glenn van de Ven 66 Benoit Famaey 55 Paolo Salucci 223344 Marco
Martorano 77 Florent Renaud 558899

(Received xxx; accepted xxx)

We present an observational study of the dark matter fraction in star-forming
disk-like galaxies up to redshift $z\sim 2.5$, selected from publicly
available integral field spectroscropic surveys, namely KMOS3D, KGES, and
KROSS. To model the $H\alpha$ kinematics of these galaxies, we employ 3D
forward-modelling, which incorporates beamsmearing and inclination
corrections, and yields rotation curves. Subsequently, these rotation curves
are corrected for gas pressure gradients, resulting in circular velocity
curves or ‘intrinsic’ rotation curves. Our final sample comprises of 263
rotationally supported galaxies with redshifts ranging from $0.6\leq z<2.5$,
stellar masses within the range
$9.0\leq\log(M_{\star}/\mathrm{M_{\odot}})<11.5$, and star formation rates
between $0.49\leq\log\left(\mathrm{SFR}\ [\mathrm{M_{\odot}\
yr^{-1}}]\right)\leq 2.5$. We estimate the dark matter fraction of these
galaxies by subtracting the baryonic contribution from intrinsic rotation
curves using dynamical models. We provide novel observational evidence,
showing that at a fixed redshift, the dark matter fraction gradually increases
with radius, indicating that the outskirts of galaxies are dark matter
dominated, similarly to local star-forming disk galaxies. This observed dark
matter fraction exhibits a decreasing trend with increasing redshift. However,
on average, the fraction within the effective radius (upto outskirts) remains
above 50%, similar to locals. Furthermore, we investigated the relationships
between the dark matter, baryon surface density, and circular velocity of
galaxies. We observe a decreasing trend in the dark matter fraction as baryon
surface densities increase, which is consistent across all stellar masses,
redshift ranges, and radii, with a scatter of 0.13 dex. On the other hand, the
correlation between the circular velocity at the outermost radius and the dark
matter fraction within this radius has a relatively low scatter (0.11 dex),
but its slope varies with stellar mass and with redshift, providing
observational evidence of the dynamical evolution of the interplay between the
baryonic and dark matter distributions with cosmic time. We observe that low
stellar mass galaxies ($\log(M_{\star}/\mathrm{M_{\odot}})\leq 10.0$) undergo
a higher degree of evolution, which may be attributed to the hierarchical
merging of galaxies.

###### Key Words.:

galaxies: kinematics and dynamics;– galaxies: late-type, disk-type and
rotation dominated; – galaxies: evolution; – galaxies: dark matter halo;–
cosmology: nature of dark matter

##  1 Introduction

According the current standard model of structure formation Dark Matter (DM),
an enigmatic form of matter devoid of any electromagnetic interaction, is
assumed to represent most of the matter in the Universe (Peebles, 1993).
Despite its significance, DM remains elusive and challenging to detect
directly. Its existence is essentially inferred from the gravitational
attraction it exerts on visible matter. In particular, flat outer galactic
rotation curves have historically served as compelling evidence for the
discrepancy between the observed dynamics and the amount of baryons, leading
to the hypothesis that galaxies are embedded within vast, extended halos of DM
that largely surpass their luminous components in mass (e.g., Rubin et al.,
1980; Bosma, 1981; Salucci, 2019, and references therein). Within the standard
cold DM framework, cosmological simulations corroborate that DM halos provide
a gravitational scaffolding that allows galaxies to form and maintain their
structures (Angulo & Hahn, 2022, and references therein), although some small-
scale challenges remain Bullock & Boylan-Kolchin (2017).

The quantity of DM within galaxies varies depending on their baryonic mass,
size, and environment, underscoring its critical influence on their overall
dynamics and evolution. Observations of nearby star-forming galaxies indeed
yield DM fraction between $40-50\%$ in their inner region within the projected
half-light radius ($R_{e}$), and between $70-90\%$ within $3R_{e}$, which
encompasses most of the galaxy and its stellar halo (Kassin et al., 2006;
Martinsson et al., 2013; Courteau & Dutton, 2015). Observations of early-type
galaxies often reveal lower DM fractions in the inner regions ($\lesssim 20\%$
within $R_{e}$, Cappellari et al. 2013), which tend to rise with increasing
stellar mass in galaxies with masses exceeding $10^{10}M_{\odot}$, but
similarly high fractions in the outskirts ($70-90\%$ within $5R_{e}$, Harris
et al. 2020). One reason may be that early-type galaxies, with their very low
gas content, have experienced extensive baryonic processes over their
lifetimes. These processes can affect both the distribution of baryons and
that of DM, as simulations have shown that baryonic processes can expel DM
over a timespan of several giga years (Pontzen & Governato, 2012). In
particular, feedback processes such as supernova explosions, stellar winds,
and active galactic nuclei stir the interstellar medium of galaxies and launch
powerful gas outflows, which induce fluctuations in mass and gravitational
potential that can in turn affect DM and diminish its fraction in the inner
region of all type of galaxies (e.g., Pontzen & Governato, 2014; Dutton et
al., 2016; Freundlich et al., 2020; Dekel et al., 2021; Li et al., 2023).
These processes can indeed dynamically heat up the DM and lead to the
formation of constant DM density cores rather than steep cusps. Therefore,
characterizing the DM fraction and its evolution with cosmic time not only
enables to better understand the influence of the DM distribution on galaxy
formation and evolution but also provides valuable insights into the physical
processes that govern galaxies and contributes to test the implementation of
feedback processes in the simulations.

In the last decade, integral field units (IFUs) in galaxy surveys have opened
up new possibilities for studying the spatially resolved kinematics of
high-$z\ $galaxies. Förster Schreiber & Wuyts (2020) notably show that by
using the resolved kinematics, it is now possible to obtain the rotation
curves (RCs) of galaxies up to $z\sim 1-2.5$. These RCs allow us to probe the
baryon and dark matter content on galactic scales, as well as their
distribution and physical properties. Genzel et al. (2017) and Lang et al.
(2017) were the first who analyzed the RCs of star-forming galaxies (SFGs) at
high-$z$ ($0.9\leq z\leq 2.4$), and found them declining; such behavior is
only seen in local massive (very high surface brightness) SFGs, while the RCs
of most normal SFGs are remarkably flat and rarely decline (e.g., Rubin et al.
1980 and Persic et al. 1996). Both studies suggest that the declining behavior
of RCs can be explained by a combination of high baryon fraction and extensive
pressure support in the inner regions. Some other high-$z$ studies of late-
type and early-type galaxies also report similar low dark matter fractions
within the effective radii (e.g., Burkert et al., 2016; Wuyts et al., 2016;
Price et al., 2016; Übler et al., 2018). Conversely, Tiley et al. (2019b)
studied the shape of rotation curves in $\approx 1500$ SFGs at $0.6\lesssim
z\lesssim 2.2$, and reported flat rotation curves similar to locals. Moreover,
Tiley et al. (2019b) reported a more than $50\%$ dark matter fraction within
$3.5\ R_{e}$, which is similar to local star-forming disk galaxies (Persic et
al., 1996; Martinsson et al., 2013; Courteau & Dutton, 2015).

In another follow-up study, Sharma et al. (2021a) studied K-band Multi-Object
Spectrograph (KMOS) Redshift One Spectroscopic Survey data, and derived the
observed kinematics using 3D-forward modelling, which is converted into
intrinsic rotation curves by correcting issues related to high-$z\
$measurements, such as beam-smearing and pressure-gradient. This provided them
with a large sample of more than 200 flat rotation curves of disk-like
galaxies at $z\sim 1$ (i.e., 6.5 Gyr look-back time). In Sharma et al.
(2021b), the authors employed the rotation curves established in Sharma et al.
(2021a) to determine the dark matter fractions using a halo model independent
technique, which estimate the invisible mass fraction needed to define
observed kinematics without invoking any dark matter halo models (e.g., NFW
Navarro et al. 1996 or Burkert Burkert 1995). In this technique, stellar mass
($M_{star}$) of the objects were estimated by fitting the spectral energy
distribution of the galaxies, and the gas (molecular and atomic) mass by means
of the scaling relations by Tacconi et al. (2018) and Lagos et al. (2011).
Assuming an exponential thin disk distribution, Sharma et al. estimated the
contribution of baryons to total mass at different scale radii (disk radius:
$0.59R_{e}$, optical radius: $1.89R_{e}$, and outer radius: $2.95R_{e}$),
where total mass is dynamical mass ($M_{dyn}\propto G^{-1}V^{2}\ R$) computed
directly from the rotation curves.111In general, the scale length (or radius)
are associated with various quantities that decrease exponentially such as the
surface brightness. The disk edge is defined as $3.2\ R_{D}$ ( $=$ $1.89\
R_{e}$), where the stellar surface luminosity $\propto\exp(-r/R_{D})$.
Therefore, scale lengths in terms of the effective radius: $R_{D}=0.59\
R_{e}$; $R_{opt}=1.89\ R_{e}$; $R_{out}=2.95\ R_{e}$.  This work showed that
the majority ($>72\%$) of star-forming galaxies in KROSS sample have dark
matter dominated ($>60\%$) outer disks ($\sim 5-10$ kpc), which agrees well
with local star-forming galaxies.

In the current study, we expand our understanding of the dark matter fraction
by incorporating a larger sample and extending the redshift range beyond that
of Sharma et al. (2021b). Specifically, we present a comprehensive
investigation using data from the KMOS3D survey, KMOS Galaxy Evolution Survey
(KGES), and previously analyzed KMOS Redshift One Spectroscopic Survey
(KROSS), covering a redshift range of $0.6<z<2.5$. Our primary objective is to
examine the redshift evolution of dark matter fractions across galactic scales
and explore questions related to the assembly history of galaxies. The
structure of this article is as follows: In Section 2, we provide an overview
of the datasets utilized in this study. Section 3 describes the employed
kinematic modeling techniques, the analysis of its outputs, and the
establishment of a robust final sample used throughout this work. Moving to
Section 4, we investigate the data for potential discrepancies, present the
results of the dark matter fraction at different galactic scales and cosmic
time, and discuss the correlations between dark matter fraction, baryonic
surface density, and circular velocity. Section 5 discusses potential caveats
associated with this study. Subsequently, in Section 6, we delve into a
detailed discussion of the main results. Finally, we summarize our findings in
Section 7. Throughout the analysis we assumed a flat $\Lambda$CDM cosmology
with $\Omega_{m,0}=0.27$, $\Omega_{\Lambda,0}=0.73$, and $H_{0}=70\
\mathrm{km\ s^{-1}}$.

Figure 1: The distributions of physical properties of KMOS3D parent and
selected sample, arranged from top-left to bottom-right: inclinations, stellar
masses, effective radii, and redshifts. The dark blue histogram represents all
$H\alpha$ datacubes with confirmed spectroscopic redshift. The blue histogram
shows $H\alpha$-detected galaxies with Flag-0 or Flag-1 in Table-6 of
Wisnioski et al. (2019). We further select galaxies with inclination angles
between $25^{\circ}\leq\theta_{i}\leq 75^{\circ}$, shown in the sky-blue
histogram. Lastly, we narrow our focus to galaxies with a signal-to-noise
ratio greater than 3, depicted by the yellow hatched histogram. This selection
criteria allow us to perform the robust kinematic modeling of these high-$z$
galaxies.

##  2 Datasets

In this study, we utilize publicly available high-$z$ galaxy surveys conducted
with the KMOS instrument (Sharples, 2014), namely the KMOS3D survey (Wisnioski
et al., 2019), the KGES (Tiley et al.,
2021)222https://astro.dur.ac.uk/KROSS/data.html. KMOS is a second-generation
spectrograph located at the Very Large Telescope (VLT), indeed offers unique
advantages for constraining dark matter fraction in high-$z$ galaxies. One key
feature of KMOS is its capability for integral field spectroscopy, enabling
efficient and rapid surveying of a large number of galaxies. Moreover, KMOS
operates in the K-band wavelength range (2.0-2.4 $\mu$m), which facilitates
the detection of rest-frame optical emission from high-$z$, such as $H\alpha$
($\lambda=0.65\ \mu$m). This wavelength coverage is furthermore of great
significance as it provides crucial information regarding the stellar
populations, star-formation rates, and kinematics of these distant galaxies.

###  2.1 KMOS3D

All objects of the KMOS3D survey are drawn from 3D-HST Survey (Brammer et al.,
2012; Skelton et al., 2014; Momcheva et al., 2016) within three extra-galactic
fields (COSMOS Scoville et al. 2007, GOODS-S van Dokkum & Brammer 2010, and
UDS Lawrence et al. 2007) covering a wide redshift range ($0.7<z<2.7$). The
KMOS3D observations in the $YJ-$, $H-$, and $K-$bands cover the $H\alpha$
emission line at redshift $0.7<z<1.1$, $1.2<z<1.8$, and $1.9<z<2.7$,
respectively. In all three filters average seeing conditions (spatial
resolution) are less than or equal to $0.55\arcsec$. We begin with analyzing
all $H\alpha$ cubes, and adopt the catalogs released in Wisnioski et al.
(2019). From the physical property catalog (Table-5, Wisnioski et al. 2019),
we mainly work with quantities such as galaxy-IDs, sky-coordinates, redshift,
magnitude, seeing, HST-axis ratios, effective radius, stellar masses, and
star-formation rates. The details, such as IDs, detection and non-detection
flags, of $H\alpha$ datacubes are given in the catalog associated with Table-6
of Wisnioski et al. (2019).

Figure 2: The distributions of physical properties of KGES parent and selected
samples, arranged from top-left to bottom-right: inclinations, stellar masses,
effective radii, and redshifts. The dark blue histogram represents all
$H\alpha$ datacubes with confirmed spectroscopic redshift. The blue histogram
shows $H\alpha$-detected galaxies Tiley et al. (2021). We further select
galaxies with inclination angles between $25^{\circ}\leq\theta_{i}\leq
75^{\circ}$, which are shown in the sky-blue histogram. Lastly, we narrow our
focus to galaxies with a signal-to-noise ratio greater than 3, depicted by the
yellow hatched histogram. This selection criteria allow us to perform the
robust kinematic modeling of these high-$z$ galaxies.

Physical properties of galaxies: Stellar masses in the KMOS3D survey
(Wisnioski et al., 2019) were derived from spectral energy distribution (SED)
modeling using Fast fitting code (Kriek et al., 2009), assuming exponentially
declining star-formation histories ($\tau>300$ Myr), solar metallicity, a
Chabrier (2003) initial mass function, and stellar population synthesis models
derived from Bruzual & Charlot (2003). The full KMOS3D sample covers the
stellar mass range $8.3\leq\log(M_{\rm star}\ \mathrm{[M_{\odot}]})\leq 11.7$.
The effective radius of galaxies are computed from high-resolution CANDELS
H-band photometry (Skelton et al., 2014), covering $-0.6\leq\log(R_{e}\
\mathrm{[kpc]})\leq 1.5$. The star-formation rate is derived using cross-
calibrated ‘ladder of SFR indicators’ (Wuyts et al., 2011). The full sample is
representative of typical star-forming galaxies (Fig. 7 of Wisnioski et al.
2019 and see Section 2.4 here below). The uncertainties in stellar masses and
star-formation rates are conservatively estimated to be 10% of the true
values. To obtain the inclination angle, we use the axis-ratio $a/b$ from
Wisnioski et al. (2019) and convert it into inclination $\theta_{i}$ according
to:

$\centering
cos^{2}\theta_{i}=\frac{(b/a)^{2}-q_{0}^{2}}{1-q^{2}_{0}}\@add@centering$ (1)

where $q_{0}$ is the intrinsic axial ratio of an edge-on galaxy (e.g. Tully &
Fisher 1977), which could in principle have a wide range of values in the
range $0.1$ – $0.65$ (e.g. Weijmans et al. 2014); here, we use the commonly
assumed value $q_{0}=0.2$, which is applicable for thick discs commonly found
at high-$z\ $(Harrison et al., 2017). The resulting inclinations cover a broad
range between $\approx 10^{\circ}-85^{\circ}$. The kinematic position angles
are adopted from Wisnioski et al. 2019, which are derived from PSF fitting
(employing Moffat-fit). The central x-y positions of sources are derived, in
this work, from H$\alpha$ cubes using sky-coordinates. We note that during the
kinematic modeling of datacubes, we remeasure the position angle and the
central x-y coordinates.

Selection criteria: The Wisnioski et al. (2019) catalog contains 627 galaxies
with spectroscopically confirmed redshift (cf. their Table 6), which we refer
to as the parent sample. From those, we select H$\alpha$-detected galaxies
(i.e., with FLAG 0,1 in the catalog) with high inclination angle
($25^{\circ}\leq\theta_{i}\leq 75^{\circ}$). This leaves us with 448 objects.
Before performing the kinematic modeling, we inspect their signal-to-
noise333The signal-to-noise ratio is defined as 10 $\times$ (signal cube/noise
cube), details of signal and noise computation can be referred from Sec.4 of
Wisnioski et al. (2019). (SN) ratio and H$\alpha$ images (same as Sharma et
al., 2021a) and accordingly divide the 448 objects of the sample in three
categories: Q1: SN$>3$ and sharp H$\alpha$ image; Q2: SN$\geq 3$ and
moderately visible source; Q3: either SN$<3$ or no appearance of the source in
the H$\alpha$ image. We discard all the Q3 galaxies, which yields a final
sample of 265 sources for kinematic modeling. The distributions of inclination
angle, stellar mass, effective radius, and redshift of the parent and selected
KMOS3D samples are shown in Figure. 1.

###  2.2 KGES

The KGES encompasses a sample of 285 galaxies located in the COSMOS, CDFS, and
UDS fields, with redshifts ranging from $1.2$ to $1.9$ Tiley et al. (2021).
The survey primarily targets the $H\alpha$, [NII]6548, and [NII]6583 emission
lines, which are redshifted to the H-band wavelength range (approximately
$1.46-1.85\mu$m). For the current study, we focus on the subset of 225 KGES
star-forming galaxies with confirmed spectroscopic redshifts, $H\alpha$
detections, and not flagged as AGNs. Additionally, we restrict our analysis to
galaxies with a K-band magnitude of K$<22.5$, similar to Gillman et al.
(2020). The median redshift of our sample is $z=1.49\pm 0.07$, and the median
seeing in the H-band observations is $0.7\arcsec$.

Physical properties of galaxies: The physical properties of the KGES sample
can be found in Tiley et al. (2019b); Gillman et al. (2020) and Tiley et al.
(2021). In particular, stellar masses were estimated by fitting the SED of
galaxies using a routine called Multi-wavelength Analysis of Galaxy Physical
Properties (MAGPHYS; da Cunha et al. 2008). The SED of each galaxy was
constructed using multi-wavelength photometry ranging from the ultraviolet to
the mid-infrared. The MAGPHYS routine compares the observed SEDs with the SEDs
from the spectral libraries of Bruzual & Charlot (2003), includes the dust
attenuation model of Charlot & Fall (2000), continuous star-formation
histories, and the initial mass function of Chabrier (2003). The resulting
stellar mass range of the $H\alpha$-detected sample is
$8.62\leq\log(M_{\text{star}}\ \mathrm{[M_{\odot}]})\leq 11.66$, as reported
in Tiley et al. (2021). The uncertainties in stellar masses are conservatively
estimated to be 10% of the true values. Star formation rates and its
uncertainties were estimated from the $H\alpha$ flux, corrected for dust
attenuation, assuming a Calzetti et al. (1994) extinction law. The estimated
star formation rates range from $0.0\leq\log(\text{SFR}\ \mathrm{[M_{\odot}\
yr^{-1}]})\leq 2.2$. Geometrical parameters such as the effective radius,
inclination, position angle, and central x-y coordinates were derived using
the GALFIT (Peng et al., 2010). For detailed calculations of all physical
quantities, we refer the reader to Gillman et al. (2020).

Selection criteria. To ensure the quality of datacubes and robustness of our
kinematic modeling, we employ the identical selection criteria that were
chosen for KMOS3D (as established in Sharma et al. 2021a). Firstly, we select
galaxies with confirmed spectroscopic redshifts detected in $H\alpha$.
Subsequently, we narrow down the sample based on two additional criteria: a)
high inclination angle ($25^{\circ}\leq\theta_{i}\leq 75^{\circ}$) and b)
$SN>3$. By applying these selection criteria, we are left with a set of 51
galaxies with sufficiently high SN ratios, which enables us to conduct
accurate kinematic modeling. Figure 2 presents the distributions of physical
properties as in Figure 1 for the KMOS3D sample.

###  2.3 KROSS: comparison sample

For comparison, we use the KROSS dataset (Stott et al., 2016), which was
previously studied in Sharma et al. (2021a, b) and Sharma et al. (2022) (and
also by Harrison et al. 2017; Johnson et al. 2018; Tiley et al. 2019a). The
KROSS targets are selected from extragalactic deep fields covered by multi-
wavelength photometric and spectroscopic data: 1) Extended Chandra Deep Field
Survey (E-CDFS: Giacconi et al. 2001; Lehmer et al. 2005), 2) Cosmic Evolution
Survey (COSMOS: Scoville et al. 2007), 3) Ultra-Deep Survey (UKIDSS: Lawrence
et al. 2007), and 4) SA22 field (Steidel et al., 1998). Some of the targets
were also selected from the CF-HiZELS survey (Sobral et al., 2015). The
targets were selected such that the $H_{\alpha}$ emission is shifted into
J-band with a median seeing of $0.7\arcsec$. The KROSS sample studied in
Sharma et al. contains 225 galaxies with redshift $0.7\leq z\leq 1.04$,
inclination range $25^{\circ}<\theta_{i}\leq 75^{\circ}$, effective radius
$0.08\leq\log(R_{e}\ \mathrm{[kpc]})\leq 0.89$, stellar mass
$9.0\leq\log(M_{*}\ \mathrm{[M_{\odot}]})<11.0$, and circular velocity
$1.45\leq\log(V_{out}\ \mathrm{[km\ s^{-1}]})\leq 2.83$, where $V_{out}$ is
calculated at $R_{out}=2.95R_{e}$.

Figure 3: Position of galaxies from star-forming main sequence (MS): we plot
offset from the MS, $\delta{\rm MS}\ (={\rm SFR}/{\rm SFR}({\rm MS};z,{M_{\rm
star}})$) as a function of the stellar mass ($M_{\star}$) and the reference
${\rm SFR}({\rm MS};z,{M_{\rm star}})$ is that of Speagle et al. (2014). The
gray shaded area represents the intrinsic scatter of 0.3 dex. The KMOS3D,
KGES, and KROSS data are denoted by red, green, and blue filled circles,
respectively. The solid lines in corresponding colors depict the average
offset from the main sequence for each dataset. All galaxies are within 1 dex
of the MS line.

###  2.4 Estimating gas masses

Observations show that typical star-forming galaxies lie on a relatively
tight, almost linear, redshift-dependent relation between their stellar mass
and star formation rate, the so-called main sequence of star formation (MS;
e.g., Noeske et al., 2007; Whitaker et al., 2012; Speagle et al., 2014). Most
stars since $z\sim 2.5$ were formed on and around this MS (e.g., Rodighiero et
al., 2011), and galaxies that constitute it, usually exhibit a rotating disk
morphology (e.g., Förster Schreiber et al., 2006; Daddi et al., 2010; Wuyts et
al., 2011). Figure 3 shows the position of the selected sources with respect
to the main sequence of typical star-forming galaxies (MS), i.e. their offset
from the main sequence

$\delta{\rm MS}={\rm SFR}/{\rm SFR}({\rm MS};z,{M_{\rm star}}),$ (2)

where ${\rm SFR}({\rm MS};z,{M_{\rm star}})$ is the analytical prescription
for the centre of the MS as a function of redshift and stellar mass proposed
in the compilation by Speagle et al. (2014), as a function of stellar mass.
The figure shows that the all sources are on or around the main sequence. This
enables us to estimate their molecular gas masses ($M_{H2}$) using the Tacconi
et al. (2018) scaling relations, which provide a parameterisation of the
molecular gas mass as a function of redshift, stellar mass, and offset from
the MS stemming from a large sample of about 1400 sources on and around the MS
in the range $z=0-4.5$ (cf. also Genzel et al. 2015 and Freundlich et al.
2019). The scatter around these molecular gas scaling relations and the
stellar mass induce a 0.3 dex uncertainty in the molecular gas mass estimates.

To calculate the atomic mass ($M_{HI}$) content of galaxies within the
redshift range $0.6\leq z\leq 1.04$, we use the HI scaling relation presented
by Chowdhury et al. (2022), which provides the first $M_{star}-M_{HI}$
relation at $z\approx 1$, encompassing 11,419 star-forming galaxies. The
relation was derived using a stacking analysis across three stellar mass bins,
each bin with a $4\sigma$ detection and an average uncertainty of $\sim 0.3$
dex. To compute the HI mass at $z>1.04$, we employ the $M_{star}-M_{HI}$
scaling relation derived from a galaxy formation model under the
$\Lambda\mathrm{CDM}$ framework (for details see, Lagos et al., 2011). This
scaling relation successfully reproduces both the HI mass functions (Zwaan et
al., 2005; Martin et al., 2010) and the ${}^{12}CO$ luminosity functions
(Boselli et al., 2002; Keres et al., 2003) at $z\approx 0$ with an uncertainty
of around $0.25$ dex, as well as follows the observations of quasars from
$z=0-6.4$ (see Fig.12 Lagos et al. 2011).

Figure 4: The distributions of physical properties of galaxies after kinematic
modelling, arranged from top-left to bottom-right: inclinations, stellar
masses, effective radii, and redshifts. These distributions are presented for
the KMOS3D (in red), KGES (in green), and KROSS (in blue) samples after
kinematic modeling. In total, we use 73 KMOS3D, 21 KGES, and 169 KROSS
galaxies to estimate the dark matter fraction at high-$z$. For detailed
information on sample selection following kinematic modeling, we refer the
reader to Section 3.2.  Figure 5: The outputs of the kinematic modeling of
selected KMOS3D (red box), KGES (green), and KROSS (blue) sources obtained
using 3DBAROLO. Columns 1-3: show the velocity map data, model, and residuals,
respectively. The gray dashed line shows the position angle, the black cross
the central x-y position, and the green line marks the plane of rotation.
Column 4: presents the major axis position-velocity diagram, where the black
shaded area with blue contour represents the data while the red contour shows
the model, and the orange squares with error bars are the best-fit line-of-
sight rotation velocity inferred by 3DBAROLO. The yellow and blue vertical
dashed lines represent the effective radius ($R_{e}$), and the optical radius
($R_{opt}=1.89\ R_{e}$), respectively. Column 5: displays the corresponding
velocity dispersion curve.

##  3 Forward modeling of the datacubes

We conduct a comprehensive reanalysis of the entire KMOS3D, KGES, and KROSS
datasets, using the 3D forward modeling approach implemented by 3DBAROLO. In
order to obtain precise kinematics, we used an optimization function in
conjunction with 3DBAROLO to more accurately constrain the essential gas
geometrical parameters, see Section 3.1. We found that galaxies with low
signal-to-noise ratios ($\sim 3-5$) pose a challenge in terms of accurate
kinematic modeling. This issue can be attributed to the effects of seeing
(PSF), which makes it difficult to constrain the kinematics of such galaxies.
Consequently, we had to discard these galaxies from the final analysis, see
Section 3.2. Our final sample consists of 73 KMOS3D, 21 KGES, and 169 KROSS
galaxies, i.e., total 263 objects. The distributions of the physical
properties of the final sample is shown in Figure 4.

###  3.1 Kinematic modeling

We model the kinematics of the galaxies in our samples using the 3DBAROLO code
(Teodoro & Fraternali, 2015). The main advantages of modeling datacubes with
3DBAROLO (hereafter 3DBarolo) are: (1) it allows us to reconstruct the
intrinsic kinematics in three spatial and three velocity components for given
initial guesses that define the kinematics and geometry of a galaxy; (2) the
3D projected modeled datacube is compared to the observed datacube in
3D-space; (3) it simultaneously incorporates the instrumental and the
observational uncertainties (e.g., spectral-smearing and beam-smearing) in
3D-space.444In the IFU based spectroscopy and imaging, the Point Spread
Function (PSF) describes how a point source (like a distant star or galaxies)
appears in an observation due to the effects of the imaging system. Similarly,
the LSF describes how a spectral line from a point source is broadened by the
spectrograph and other instrumental effects. That is, beamsmearing is
associated to poor PSF, while spectral smearing is due to poor LSF. For
details, we refer the reader to Teodoro & Fraternali (2015) and Di Teodoro et
al. (2016). This 3-fold approach of deriving kinematics is designed to
overcome the observational and instrumental effects and hence allows us to
stay close to the realistic conditions of the galaxy. Therefore, it gives us
somewhat improved results than the 2D-approach555The 2D-kinematic modeling is
a technique that uses datacubes to create detailed velocity-maps and then
rotation curve is derived along the projected major axis of the galaxy.,
specifically, in the case of small angular sizes and moderate SN of high-$z$
galaxies (see Di Teodoro et al. 2016). Basic assumptions under 3DBarolo, its
basic requirement, and limitations are detailed in Sharma et al. (2021a,
section-3.1), and briefly mentioned below.

The kinematic modeling with 3DBarolo requires three geometrical parameters–
namely the galaxy’s central position ($x_{c},y_{c}$), the inclination angle
($\theta_{i}$), and the position angle (PA), and three kinematic parameters–
namely the redshift ($z$), the rotation velocity ($V_{rot}$), and the velocity
dispersion of the ionized gas ($V_{disp}$). In our modelling, we fix the
geometrical parameters and redshift, while kinematic parameters are free to
determine. 3DBarolo comes with several useful features particularly useful for
high-$z$ low S/N data666Cf. its documentation,
https://bbarolo.readthedocs.io/en/latest/. We use 3DFIT TASK for performing
the kinematic modelling. 3DBarolo produces mock observations given the input
parameters in the 3D observational space ($x,y,\lambda$), where ($x,y$) stands
for the spatial axes and $\lambda$ is the spectral axis coordinate, resulting
in a datacube: $f(x,y,\lambda)$. These models are fitted ring-by-ring to the
observed datacube in the same 3D space, accounting for beam smearing. A
successful run of 3DBarolo delivers the beam smearing corrected velocity (or
moment) maps, the stellar surface brightness profile, the rotation curve (RC),
and the dispersion curve (DC) along with the kinematic models. Notice, the RC
(or position velocity-diagrams) are not derived from the velocity maps, but
instead calculated directly from the datacubes by minimizing the difference in
$V_{LOS}=V_{rot}\sin\theta_{i}$ of model and data in each ring.

When we employed 3DBarolo to model the galaxies, we discovered that the
geometrical parameters ($x_{c},y_{c}$, and PA) were inadequate for accurately
representing their kinematics. The reason being that gas kinematics differ
from stellar kinematics due to their different morphologies (i.e., geometry).
In particular, as we go higher in the redshift gas morphological parameters:
PA, $x_{c}$, and $y_{c}$, starts differentiating from their photometric
measurements, i.e., stellar morphology (as reported in Wisnioski et al., 2015;
Harrison et al., 2017; Sharma et al., 2021a). Moreover, from our previous work
with 3DBarolo (Sharma et al., 2021a), we also learned that 3DFIT TASK is
incapable to constrain simultaneously multiple (more than 2) parameters.
Therefore, in this work, we estimated the gas geometrical parameters using an
optimization function that runs atop 3DBarolo, namely minimizing the following
loss function:

$\ell=\sqrt{\overline{(X_{D}-X_{M})^{2}}}\times\frac{N[D]}{N[D:\neq
0,\neq\pm\infty]}\times\frac{N[M]}{N[M:\neq 0,\neq\pm\infty]}$ (3)

where $X_{D/M}$ is array of data/model, N represents the length of the array,
and D & M stands for data & model, respectively. In equation, the first term
is root mean square error, and last two terms are the weights from data and
model. In denominator, $N[D:\neq 0,\neq\pm\infty]$ gives the length of
datacube which includes only non-zeros and finite elements, same for model
$N[M:\neq 0,\neq\pm\infty]$. That is, weights are higher if datacube contains
more zero and infinite elements; hence, high value of loss function ($\ell$).
We use a Nelder-Mead minimization method777Nelder-Mead is a gradient free
optimization method that finds the minimum of function by iteratively updating
the vertices of polytope in n-dimensional space. It is particularly used for
optimization problem with non-liner or multi-model objective functions, which
is the case of galaxy kinematics., which is available in the scipy.optimize
library. This optimization function enables to fit multiple parameters with
3DBarolo. We remark that, we ran the 3DBarolo on each object as we did in
Sharma et al. (2021a), i.e., free parameters in 3DBarolo are only
$v_{rot},v_{disp}$; while the extra parameters are constrained by optimizer.
The details of logical flow of optimization function with 3DBarolo is
described in Appendix A.

###  3.2 Inspection of kinematic modeling outputs

We modeled the kinematics of 541 star-forming galaxies, those include 265 from
KMOS3D, 51 from KGES, and 225 from KROSS. For quality assessment and
assurance, we inspected the outputs of 3DBarolo+optimization for each
individual galaxy. Firstly, we scrutinized the optimization log of all KMOS3D
and KGES galaxies. In KMOS3D, we notice that, 57 of them experienced
optimization failure due to low SN. In addition, 3DBarolo could not perform
the modeling for 52 galaxies owing to their large PSF. Furthermore, we
observed that 78 galaxies had a maximum radius comparable to the PSF, thereby
allowing Barolo to form only two rings, the first of which being unreliable
(see Di Teodoro et al. 2016 and Sharma et al. (2021a)). Consequently, rotation
curves obtained from these galaxies can not be used in dark matter fraction
studies. Lastly, we inspected the velocity maps and high-resolution
photometric images, and we notice 5 galaxies have disturbed kinematics due to
nearby neighbors, therefore discarded. After kinematic modeling, in total, we
had to exclude 192 objects. It is not surprising for us to lose a lot of
high-$z$ data, as these observations are often noisy and angular size of
objects is very small. This is consistent with the findings of previous
studies on KMOS3D data, which were often conducted on a relatively small sub-
sample (see, Genzel et al., 2017, 2020; Nestor Shachar et al., 2022, and the
references therein). The final KMOS3D sample contains 73 galaxies, which is a
large enough sample to perform a statistical study.

Within the KGES dataset, three galaxies encountered optimization failures,
while 10 galaxies exhibited a larger PSF compared to the actual galaxy size.
Additionally, 16 galaxies displayed a maximum radius that was comparable to
the PSF, thereby allowing 3DBarolo to form only two rings, the first of which
being unreliable. Consequently, these galaxies were deemed unsuitable and were
excluded from further analysis. As a result, our final KGES sample contains
only 21 galaxies. Furthermore, for the details of comparison sample, i.e.,
KROSS objects, see Appendix A (Para.-4). In the end, we retain 263 galaxies
(73 from KMOS, 21 from KGES and 169 from KROSS).

Figure 6: The impact of the pressure gradient on the circular velocity of
galaxies computed at $R_{out}$. The KMOS3D, KGES, and KROSS samples are
denoted by red circles, green cross, and blue stars, respectively. The
interior of each data point is color-coded for rotation-to-dispersion ratio
computed before pressure support correction ($V^{int}_{{Rout}}/\sigma$). The
black dashed line represents the one-to-one relation between the two
measurements. We observe that galaxies with higher velocity dispersion require
large pressure corrections.

###  3.3 Kinematic modeling results

Figure 5 shows a few examples of the kinematic modeling results. In this
figure, the first three columns in each row show the moment-1888The first
moment (or moment-1) is the intensity-weighted velocity map, which shows the
mean velocity of the emission cube at each point., data, model, and residuals,
respectively, from left to right. As we can see, the data closely matches the
model, and residuals lie around zero, which was precisely the aim of using
both 3DBarolo and the optimization function. The fourth column displays the
major axis position-velocity diagram (hereafter PV-diagram, which is the line-
of-sight (LOS) rotation velocity of the galaxy at each spatial-bin.999Spatial-
bin is based on the size of the spherical-ring in which line-of-sight velocity
is computed. The red contour represents the model, and the black shaded area
with blue contour represents the data. The orange squares with error bars
indicate the best-fit LOS rotation velocity. The yellow and blue vertical
dashed lines represent the effective radius ($R_{e}$) and optical radius
($R_{opt}=1.89\ R_{e}$) of the galaxy, respectively. The last column presents
the velocity dispersion curve.

To estimate the DM fraction, we first correct the rotation curves for pressure
support, as the turbulent interstellar medium in high-$z$ galaxies makes them
pressure-supported systems (Burkert et al., 2010; Übler et al., 2019). In
Sharma et al. (2021a), it was demonstrated that high-z galaxies exhibit a non-
uniform and non-isotropic velocity dispersion, inducing a pressure gradient.
This pressure gradient significantly hampers the motion of gas, leading to a
decrease in the rotation velocity of the gas in the inner region of galaxies
by $\sim 50\%$ of its original value. In some cases, it also affects the outer
rotation curves, causing them to decline.

To address this issue, Sharma et al. (2021a) introduced a method known as the
‘Pressure Gradient Correction’ (PGC), which effectively corrects for the
pressure support. For a detailed explanation of this method, we refer the
reader to (Sharma et al., 2021a, Sec. 3.2). In this study, we applied the PGC
method to all datasets and investigated the impact of velocity dispersion on
the circular velocity of galaxies. In Figure 6, we present the intrinsic
velocities ($V^{int}$) and pressure-corrected circular velocities ($V^{PGC}$)
within $R_{out}$, color coded for intrinsic rotation-to-dispersion ratio
($V^{int}_{Rout}/\sigma$), where $V^{int}$ is rotation velocity without
pressure support corrections. As shown, systems primarily supported by
rotation ($V^{int}_{Rout}/\sigma>1$) exhibit minimal or no pressure
correction, whereas dispersion-dominated systems induce significant
correction. Additionally, there are substantial pressure corrections at the
lower end of the velocity range ($\log(V_{R_{out}}[km/s])<2.3$), that
gradually decrease towards higher velocities and approaches to nearly zero. We
note the readers that, prior to the implementation of PGC, there were only 9
dispersion dominated galaxies (3 KMOS3D, 1 KGES, and 5 KROSS). However, after
applying PGC, none of these galaxies have $V^{int}_{Rout}/\sigma<1$, as
depicted in Figure 16. Therefore, we do not exclude these galaxies from our
analysis. Hence, the full sample is a good representative of rotation
supported system.

We finally examined the relationship between the size of the point spread
function (PSF) and the effective radius and found that 59% (43) of KMOS3D, 90%
(19) of KGES, and 88% (148) of KROSS galaxies possess a PSF larger than their
effective radius. This suggests that the majority of the sample cannot be used
for studying the dark matter fraction within $R_{e}$, as doing so would result
in highly uncertain outcomes. Below, we only use the 53 galaxies with a large
enough effective radius compared to the PSF to characterize the DM fraction
within that radius. Additionally, 16% (12) of KMOS3D, 19% (4) of KGES, and 42%
(71) of KROSS galaxies have a PSF larger than their optical radius. However,
1.4% (1) of KMOS3D, 0.0% (0) of KGES, and 12% (20) of KROSS galaxies have a
PSF larger than their outer radius. Therefore, the most reliable measurement
of the dark matter fraction is obtained from the outer radius ($R_{out}$).
Consequently, in this work we mostly focus on interpreting the results that
are computed within (or at) $R_{out}$.

##  4 Results

Our primary goal is to ascertain the dark matter fraction across various
galactic scales. However, recent studies at high-$z\ $suggest that galaxies
contain baryon-dominant inner regions ($<R_{e}$), where dark matter
constitutes less than 20% of the total mass (Genzel et al., 2020; Nestor
Shachar et al., 2022). Adding to this perspective, a study by Lelli et al.
(2023) found that the kinematics of two main-sequence galaxies, from the
cosmic dawn, could be entirely explained by baryonic dynamical mass-models,
eliminating the need for dark matter. These findings amplify the long-standing
disk-halo degeneracy issue, a challenge that persists even in local galaxies
(van Albada et al., 1985), and remains unresolved, particularly when mass-to-
light ratios are not entirely reliable (Sofue & Rubin, 2001; Bullock & Boylan-
Kolchin, 2017). Therefore, in this work, we adopt a conservative approach,
beginning with the analysis of rotation curves that contemplate the conditions
of maximal baryonic disk. In particular, we mass-model the rotation curves
under two scenarios that consider the maximum contribution from: (1) stellar
disk, and (2) baryonic (stars+gas) disk. The modeling details are compiled in
Appendix B.

In Figure 7, we present a comparison between the fiducial and the mass-
modelled stellar and baryonic masses in the case of maximum stellar disk and
maximum baryonic disk scenarios (upper and lower panels, respectively). When
considering the maximum stellar disk, we observe that fitting the rotation
curves requires stellar masses that are twice the values obtained from
photometric measurements for the majority of the sample. This finding strongly
suggests the requirement of an additional gas component at high-$z\ $, which
we incorporate in the maximum baryonic disk scenario. However, as illustrated
in the bottom panel of Figure 7, even after including the gas component, 50%
of the sample still demands baryonic masses that are factor of two higher than
their fiducial values. Such values are unrealistic to obtain within the range
of observed uncertainties. That is, with a large sample of 263 galaxies, we
are unable to rule out the presence of dark matter halos at high-$z\ $.

Moreover, recent observations of high-$z\ $galaxies suggest that longer
integration time is crucial for accurately mapping the complete stellar and
gas kinematics. In particular, a work by Puglisi et al. (2023) on KMOS Ultra-
deep Rotation Velocity Survey (KURVS) has demonstrated that deep observations
enhance both the amplitude and radial extent of the rotation curves.
Consequently, deep observations of our current sample will further demand the
inclusion of an additional halo component, as anticipated while analyzing the
mass-modeled rotation curves (see Figure 18). Henceforth, we proceed the
investigation of the rotation curves in the presence of a dark matter halo
component.

Figure 7: Comparison of photometric (PHOT) stellar and baryonic (stars+gas)
masses with mass-modeled (MM) stellar and baryonic masses in the case of the
maximum stellar disk and maximum stellar+gas disk scenarios, upper and lower
panels, respectively. The KMOS3D, KGES, and KROSS datasets are depicted in
red, green, and blue colors, respectively. The dashed black line indicates
when mass-modeled masses are equivalent to their fiducial values (i.e.,
photometric masses), suggesting that objects lying above this line require the
presence of a dark matter halo. Figure 8: Dark matter fraction within $R_{e}$
as a function of stellar masses for galaxies having $R_{e}>R_{\mathrm{psf}}$.
The KMOS3D, KGES, and KROSS datasets are depicted in red, green, and blue
colors, respectively. For comparison, we include two local massive disk
galaxies: the Milky Way (MW) and Andromeda (M31), represented by the blue and
the orange star, respectively. Additionally, we compare these results with
previous studies at high redshift, such as Genzel et al. (2017, yellow
circles), Genzel et al. (2020, yellow hexagons), Nestor Shachar et al. (2022,
yellow dots), and Bouché et al. (2022, yellow star). The black square
represents the two galaxies that are in-common with Genzel et al. (2020) drawn
from this work, while pink hexagons shows the measurements of same objects
from Genzel et al. (2020). The black and yellow horizontal dashed lines
represent the regimes of 100% dark matter dominance and baryon dominance,
respectively. The gray shaded area indicates the forbidden region, where
galaxies with $M_{\mathrm{dyn}}<M_{\mathrm{bar}}$ are located.  Figure 9: Dark
matter fraction within $R_{out}$. The KMOS3D, KGES, and KROSS datasets are
represented in red, green, and blue colors, respectively. For comparison, we
include two local massive disk galaxies: the Milky Way (MW) and Andromeda
(M31), represented by the blue and the orange star, respectively. The pink
shaded area represents local star-forming disks (Persic et al., 1996). The
solid black dashed line represents the 100% dark matter regime, while the
yellow dashed line represents the baryon-dominated regime. The gray shaded
area indicates the forbidden region, where galaxies with $M_{dyn}<M_{bar}$ are
located. These color codes are same through out the text.

The total dark matter mass can be calculated by subtracting the baryonic mass
contribution from observed rotation curve, which gives the estimate of total
dynamical mass ($V^{2}(R)=\sqrt{GM_{dyn}/(\kappa R)}$ where $\kappa$ is a
geometric correction to spherical symmetry) as a function of radius, as
explained in Appendix B. That is, dark matter halo modeling is not necessary
to study the amount of dark matter. In order to estimate the dark matter
fraction of our datasets, we use Equation 9, which is a halo model independent
approach (previously postulated in Sharma et al. 2021b), yet allow us to
estimate the dark matter at different galactic scales. We begin by inspecting
dark matter fraction within $R_{e}$, which is a rather controversial issue for
high-$z$ galaxies. For example, Genzel et al. (2020) and Nestor Shachar et al.
(2022) reports dark matter-deficiency within $R_{e}$ at high-$z\ $, on the
other hand, Sharma et al. (2021b) and Bouché et al. (2022) reported a similar
amount ($>45\%$) of dark matter fraction as seen in local star-forming disk
galaxies.

###  4.1 Dark matter fraction within $R_{e}$

Here, we present the dark matter fraction only for the galaxies that have
$R_{e}>R_{psf}$. In total, we have only 20% (53) galaxies that abide this
criteria, shown in Figure 8. We observe that five galaxies fall in the
forbidden region, indicating either inaccurate photometric stellar mass
estimates or incomplete sampling of stellar and gas motion in the observed
rotation curves. Apart from these five forbidden galaxies, only five
additional galaxies exhibit $f_{\mathrm{DM}}(<R_{e})<50\%$, while the
remaining galaxies in the sample are dark matter dominated within $R_{e}$.
Notably, the dark matter fraction of very massive ($>10^{10.5}\
\mathrm{M_{\odot}}$) galaxies in our sample surpasses that of local massive
galaxies, such as the Milky Way (Petač, 2020, and ref. therein) and Andromeda
galaxies (Tamm et al., 2012, and ref. therein), represented by blue and orange
stars in Figure 8, respectively.

We compared our findings with those of previous studies, such as Genzel et al.
(2017, 2020); Nestor Shachar et al. (2022), and Bouché et al. (2022), as shown
by the different shapes of yellow markers in Figure 8. Our results are in full
agreement with Bouché et al. (2022), and a few galaxies of Genzel et al.
(2020) and Nestor Shachar et al. (2022). However, we observe that the majority
of Nestor Shachar et al. (2022), Genzel et al. (2020), and Genzel et al.
(2017) are baryon-dominated. Specifically, the very massive galaxies studied
in Genzel et al. (2020) and Nestor Shachar et al. (2022) are dark matter
deficient, which contrasts with our findings, where galaxies in the same mass
range have $\sim$50%-80% dark matter.

Next, we cross-matched KMOS3D sample and the galaxies studied in Genzel et al.
(2020), resulting in the identification of only two overlapping systems:
GS4_05881 ($z=0.99$) and GS4_43501 ($z=1.61$). These systems are represented
by black squares in Figure 8. According to Genzel et al. (2020), the reported
dark matter fractions within $R_{e}$ for these two systems are
$f_{{}_{DM}}(<R_{e})=0.64$ and $0.19$, respectively, shown by pink hexagons.
While, our estimates are about $0.98$ and $0.75$, respectively. These values
are $\approx 1.5$ and 4 times higher, respectively, than the estimates
reported in Genzel et al. (2020). This difference is most-likely attributed to
the distinct kinematic modeling and pressure support corrections implemented
in this work, quantified and reported previously in Sharma et al. (2021a, b).
Moreover, it is noticeable that the stellar masses of these objects reported
in Genzel et al. (2020) are marginally higher than those in the present study.
This discrepancy arises because we employ photometric stellar masses, whereas
Genzel et al. (2020) utilizes dynamically mass-modeled (best-fit) stellar
masses.

Furthermore, we had the opportunity to refine the measurement of $R_{e}$ for a
subset of the KROSS sample through the latest observations from the James Webb
Space Telescope (JWST). This particular subset of galaxies is situated within
the COSMOS field (Skelton et al., 2014), which has recently undergone
observation by the COSMOS-WEB team (Casey et al., 2023). Estimates of $R_{e}$
for these galaxies were derived using GALFITM (Häußler et al., 2013) setup
equal to that presented in Martorano et al. (2023). As discussed in Appendix D
and illustrated in Figure 19, we observe that the new $R_{e}$ estimates shift
the galaxies of low dark matter fractions ($<20\%$) towards higher values. In
this specific sub-sample, none of the galaxies have a dark matter fraction
below $20\%$. Consequently, we suggest that a low dark matter fraction within
$R_{e}$ is very unlikely for disk-like galaxies at high-$z\ $($0.65\leq z\leq
2.2$).

Figure 10: Dark matter fraction across the cosmic time ($0.0\leq z<2.5$). Left
panel: dark matter fraction within $R_{out}$ as a function of redshift. The
color codes of datasets are same as Figure 9, and given in the legend of the
plot. The gray stars represents the binned data points with 1$\sigma$
uncertainty shown by error-bars. Right panel: the averaged dark matter
fraction within $R_{e}$ and $R_{out}$ as a function of redshift, represented
by yellow and blue lines, respectively. The bootstrap-errors are indicated by
shaded areas color-coded with the same color as the averaged line. Dark matter
fraction of local star-forming galaxies within $R_{out}$ are from Persic et
al. (1996) and within $R_{e}$ are taken from disk galaxy survey Courteau &
Dutton (2015). The uncertainties on these measurements are, on average, $\pm
0.25$. The color code of these estimates are same as high-$z$ galaxies. The
dark-yellow dotted-dashed lines shows the dark matter fraction within $R_{e}$
from Nestor Shachar et al. (2022).

###  4.2 Dark matter fraction within $R_{out}$

In Figure 9, we plot the dark matter fraction within $R_{out}$ ($\sim 3\times
R_{e}$), as a function of stellar mass. Firstly, we observe that the majority
of the sample exhibits dark matter-dominated outer disks, with 16% of objects
showing $0.2\leq f_{\mathrm{{}_{DM}}}(<R_{out})<0.5$, and only 6% of objects
having $f_{\mathrm{{}_{DM}}}(<R_{out})<0.2$ (including objects in the
forbidden region). Secondly, we notice a slightly decreasing trend of
$f_{{}_{DM}}$ as a function of stellar mass, except for a few outliers. These
outliers are six massive galaxies in the KMOS3D sample that display prominent
disks in HST images and exhibit rising-flat rotation curves, as shown in
Figure 23.101010The following massive galaxies have very high dark matter
fractions: COS4_23890 ($z=0.85$, $M_{*}=10^{9.83}$, $f_{{}_{DM}}=0.90$),
U4_19708 ($z=1.66$, $M_{*}=10^{11.26}$, $f_{{}_{DM}}=0.83$), U4_31577
($z=1.52$, $M_{*}=10^{11.07}$, $f_{{}_{DM}}=0.85$), and U4_26875 ($z=1.36$,
$M_{*}=10^{11.18}$, $f_{{}_{DM}}=0.86$), U4_28358 ($z=1.42$,
$M_{*}=10^{10.98}$, $f_{{}_{DM}}=0.84$), U4_21533 ($z=0.92$,
$M_{*}=10^{11.11}$, $f_{{}_{DM}}=0.71$). Notably, their dark matter halos are
more massive than those of the most massive systems in the present-day
Universe, such as the Milky Way (MW) and Andromeda (M31), shown by blue and
orange stars, respectively, in Figure 9.

In Appendix D, Figure 20 shows the results of dark matter within $R_{opt}$
($\sim 2\times R_{e}$). We notice that dark matter fraction is on average 10%
less within $R_{opt}$ than $R_{out}$, suggesting dark matter dominates the
outer-disks at high-$z$, which is very similar to local disk galaxies (Persic
et al., 1996). Finally, we investigate galaxies with
$f_{\mathrm{{}_{DM}}}(<R{out})<20\%$ in relation to their PSF and SN, but no
dependencies are found. To explore further, we estimated the stellar masses of
these objects using the S$\acute{e}$rsic profile. It is important to note that
the stellar mass derived from the Freeman disk assumption is, on average, 1.02
times higher than when assuming the S$\acute{e}$rsic profile. Consequently,
the dark matter fraction derived using the S$\acute{e}$rsic profile is only 2%
higher (see Appendix D and Figure 21). Therefore, altering the assumed stellar
mass profile/distribution does not significantly improve the estimated dark
matter fraction for these objects that fall into the so-called forbidden
region or show low $f_{\mathrm{{}_{DM}}}(<20\%)$.

###  4.3 Dark matter fraction across cosmic time

To gain a deeper understanding of the evolution of dark matter with cosmic
time, we initially divided the dark matter fraction,
$f_{\mathrm{{}_{DM}}}(<R{out})$, into three redshift bins: $0.5<z\leq 1.1$,
$1.1<z\leq 1.8$, and $1.8<z\leq 2.5$. We employed weighted mean statistics to
determine the binned values, and the errors were estimated using a bootstrap
analysis. The results are presented in Figure 10 (left panel), where the
binned data points are denoted by large gray stars connected by a solid black
line. Upon inspecting the figure, it becomes evident that the dark matter
fraction exhibits a decreasing trend between redshift range of $0\leq z\leq
1.8$. On the other hand, above redshift 1.8, there appears to be a potential
increase in the dark matter fraction. However, the latter is subject to
debate, and in order to confidently establish this trend, additional data is
required in the high-redshift bin ($z>1.8$). Consequently, we refrain from
further discussing the higher-redshift bin ($z>1.8$) in our analysis.

In the right panel of Figure 10, we present the binned dark matter fraction
within $R_{e}$ and $R_{out}$, as a function of redshift. To provide context
and comparisons, we have included measurements from local star-forming
galaxies. In particular, we used the disk galaxies sample from Persic et al.
(1996), their results align closely with those obtained from the disk galaxy
survey conducted by Courteau & Dutton (2015). As shown in figure, the dark
matter fraction exhibits a decreasing trend between $0\leq z\leq 1.8$ at both
radii. Our findings regarding $f_{{}_{\mathrm{{DM}}}}(<R_{e})$ stand in
contrast to those of Nestor Shachar et al. (2022), who report a rapid decrease
in $f_{{}_{\mathrm{{DM}}}}(<R_{e})$ as a function of redshift. Statistically,
the dark matter fraction within the effective radius does not go below 50%.
However, it is worth noting that in order to further refine our understanding
on dark matter content within $R_{e}$, more resolved (high-quality)
observations are required.

Figure 11: Scaling relations of dark matter, $f_{{}_{DM}}-\Sigma_{bar}$ and
$f_{\mathrm{{}_{DM}}}-V_{c}$, shown in the left and right panels,
respectively. The blue dashed line represents the best fit to the data in both
relations, with the associated intrinsic scatter ($\zeta_{int}$) indicated on
the plot. In the $f_{\mathrm{{}_{DM}}}-\Sigma_{bar}$ relation, a comparison
fit from Genzel et al. (2020) is shown as a red dashed line. The color codes
used for data in both panels are consistent with Figure 9.

Figure 12: Exploration of the $f_{{}_{DM}}-\Sigma_{bar}$ relation across
different stellar masses and redshift intervals. Upper panel:
$f_{{}_{DM}}-\Sigma_{bar}$ relation depicted for the low-stellar mass bin
($M_{\text{star}}\leq 10^{10}\ M_{\odot}$) in yellow and the high-stellar mass
bin ($M_{\text{star}}>10^{10}\ M_{\odot}$) in orange. The blue dashed line
represents the best fit to the data given by Equation 4. Additionally, a
comparison fit from Genzel et al. (2020) is shown in red dashed line. Lower
panel: $f_{{}_{DM}}-\Sigma_{bar}$ relation separated for low- and high-stellar
masses, left and right panels respectively. In both panels, low-$z$ ($z\leq
1$) and high-$z$ ($z>1$) galaxies are shown in blue and red, respectively, and
corresponding best fits are shown in dark blue and red dashed-lines. In all
panels, median error on the data points is indicated by gray crosses, and
filled crosses & circles represent the measurements taken within $R_{opt}$ and
$R_{out}$, respectively. We note the reader that these relations are plotted
exclusively for objects with $f_{{}_{\rm DM}}>20\%$, while the full sample is
shown in Figure 11.

###  4.4 Dark matter scaling relations

Here, we present the dark matter fraction correlation with the baryon surface
density– $\Sigma_{\text{bar}}(<R_{\text{out}})$ and the circular velocity–
$V_{\text{Rout}}$ of galaxies, as shown in the left and right panels of Figure
11, respectively. We fit both correlations using an exponential power law,
represented mathematically as follows:

$\displaystyle f_{{}_{DM}}(<R_{\text{out}})$ $\displaystyle=$ $\displaystyle
1-10^{-0.34+0.70(\log(\Sigma_{\text{bar}})-8.5)}$ (4) $\displaystyle
f_{{}_{DM}}(<R_{\text{out}})$ $\displaystyle=$ $\displaystyle
1-10^{-1.36(\log(V_{c})-1.71)}$ (5)

These relations are fitted using the ”scipy curve fit” module, by minimizing
the intrinsic scatter.

As shown in Figure 11, the dark matter fraction displays a negative
correlation with baryon surface density, which is expected and observed at
lower redshifts (see McGaugh 2010). While, it exhibits a positive correlation
with circular velocity. Together these relations imply that, although the
dynamics of galaxies are dominated by dark matter, the baryons still play an
important role in hampering the presence of dark matter, i.e., the
evolutionary stages of baryonic matter most likely seem to strongly impact the
distribution of dark matter within galaxies. This has long been known at low
redshift, but it is highly interesting that the trend seems to continue at
higher redshift, too.

To thoroughly examine the scenario postulated above, we divided our full
sample into two stellar mass bins: a low-mass bin ($M_{\text{star}}\leq
10^{10}\ M_{\odot}$) and a high-mass bin ($M_{\text{star}}>10^{10}\
M_{\odot}$). We plotted the
$\Sigma_{\text{bar}}(<R_{\text{out}})-f_{{}_{\text{DM}}}(<R_{\text{out}})$
correlations in the upper panel of Figure 12. We note the reader that this
relation is plotted exclusively for objects with $f_{{}_{\rm DM}}>20\%$. The
low stellar mass objects are represented by bright yellow, while the high mass
objects are depicted in orange. We observed that the
$\Sigma_{\text{bar}}(<R_{\text{out}})-f_{{}_{\text{DM}}}(<R_{\text{out}})$
relation, as given in Equation 4 (and Genzel et al., 2020), remains the same
for both mass ranges, indicating no change in the slope. Moreover, this
relation remains consistent when observed for $f_{{}_{\rm DM}}(<R_{opt})$,
represented by cross marks, so it is valid at different radii within galaxies.

Conversely, when we plot $V_{\text{Rout}}-f_{{}_{\text{DM}}}$ (at both radial
scales: $R_{\text{opt}}$ and $R_{\text{out}}$) for different mass bins, we
observed a distinct offset between low and high stellar mass galaxies, as
shown in the top panel of Figure 13. In particular, galaxies with high stellar
mass are fast-rotating systems with a relatively low dark matter fraction at
the outer radius for a given rotational velocity, while the opposite trend is
observed for low stellar mass systems. This likely suggests an evolution in
the distribution of dark matter due to baryonic processes that take place in
massive galaxies. Next, we segregated low and high stellar mass galaxies into
low and high redshift ($z\leq 1$ and $z>1$, respectively) as shown in the
bottom panel of Figure 13. We observed very clearly distinct behaviors within
the low mass galaxy population. Specifically, at $z>1$ in this population,
galaxies exhibited lower dark matter fractions at a given rotational velocity
compared to the same mass range at $z\leq 1$. In contrast, high mass systems
demonstrate a more similar behavior for the two redshift bins, with only a
slight variation in circular velocities and dark matter fractions.111111We
note the reader that correlations in Figure 12 & 13 are plotted exclusively
for objects with $f_{{}_{\rm DM}}(<R_{opt})>20\%$.

We also segregated low and high stellar mass galaxies into low and high
redshift ($z\leq 1$ and $z>1$, respectively) for the baryon surface density
relation, as shown in the bottom panel of Figure 12. We notice that in the
high stellar mass bin, both low and high-$z$ galaxies follows a similar trend
with the same shape as given by Equation 4. The same is true for the lower-$z$
galaxies of the low stellar mass bin. However, we noted that the high-$z$
galaxies of the low stellar mass bin display a sharp cutoff at surface
densities lower than $log(\Sigma_{\text{bar}}(<R_{\text{out}}))\approx 8.0$.
This is most likely hinting an observational bias, i.e., low baryon density
galaxies seem to be missing at $z>1$. This causes the fit to not follow the
Equation 4, but the discrepancy is not as clear visually as for the
$V_{\text{Rout}}-f_{{}_{\text{DM}}}(<R_{\text{out}})$ relation where the
redshift evolution is obvious.

We finally note that the correlation between
$\Sigma_{\text{bar}}(<R_{\text{out}})$ \- $f_{{}_{DM}}(<R_{\text{out}})$
exhibits, in total, a higher intrinsic scatter ($\zeta_{\text{int}}=0.13$ dex)
compared to the correlation between $V_{\text{Rout}}$ \-
$f_{{}_{DM}}(<R_{\text{out}})$ ($\zeta_{\text{int}}=0.11$ dex), but some
systematic uncertainties (see next Section) can artifically increase those
scatters : hence, these should not be taken as definitive before having higher
quality data. The trends, on the other hand, are clear and solid.

Figure 13: Exploration of the $f_{{}_{DM}}-V_{c}$ relation across different
stellar masses and redshift intervals. Upper panel: $f_{{}_{DM}}-V_{c}$
relation for the low-stellar mass bin ($M_{\text{star}}\leq 10^{10}\
M_{\odot}$) in yellow and the high-stellar mass bin ($M_{\text{star}}>10^{10}\
M_{\odot}$) in orange. The best-fit to original $f_{{}_{DM}}-V_{c}$ relation
is shown blue dashed line given by Equation 5. The Equations corresponds to
best fit of low- and high-stellar mass bins is written on the plot, shown by
brown and orange dashed lines, respectively. The intrinsic scatter
($\zeta_{int}$) associated to fits are printed on plot using same color code.
Lower panel: $f_{{}_{DM}}-V_{c}$ relation separated for low- and high-stellar
masses, left and right panel respectively. In both panels, low-$z$ ($z\leq 1$)
and high-$z$ ($z>1$) galaxies are shown in blue and red, respectively, and
corresponding best-fits are shown by dashed blue and red lines. In all panels,
filled crosses and circles represent the measurements taken within $R_{opt}$
and $R_{out}$, respectively, and median error on the data points is indicated
by gray crosses. Low surface brightness and low velocity low mass galaxies at
high redshift might be missed from the sample. It is, on the other hand,
however clear that DM fractions are lower at fixed rotational velocity at
higher redshift, i.e. the plots are not biased on the vertical axis. We note
the reader that these relations are plotted exclusively for objects with
$f_{{}_{\rm DM}}>20\%$ for the clarity on the matter, while the full sample is
shown in Figure 11.

##  5 Caveats

We estimate the dark matter fraction using a halo model-independent approach,
assuming an exponential thin disk distribution for stars and gas. Our analysis
provides reliable measurements of the dark matter fraction within different
radii: $R_{e}$, $R_{opt}$, and $R_{out}$. However, it is important to note
that only 20% (53 objects) of the total sample (263 objects) satisfy the
criterion $R_{e}>R_{psf}$. Nevertheless, we report that galaxies are not dark
matter deficient at high-$z$, especially, within $1-3R_{e}$. Although we
report a high dark matter fraction within $R_{e}$ for the high-$z\ $galaxies
in a larger sample compared to previous studies, we advise the reader to
interpret these measurements with caution. For instance, the JWST photometry
will have the capability to resolve the inner stellar disk and bulge
components of galaxies, for example, see Section D and Figure 19. Moreover,
deep spectroscopic observations will further improve the shape of the rotation
curves, as reported in (Puglisi et al., 2023, Fig. 2). Both the former and the
latter will impose even tighter constraints on the dark matter distribution in
the inner regions of galaxies, which may differ from the results presented
here. In particular, the geometry of the disk, especially the gaseous disk,
could be different (Renaud et al., 2021), whilst many systematic uncertainties
in the actual stellar mass and gas mass could change the results a bit,
although we have shown in Sect. 4.1 that it would take a very serious mismatch
to erase the dark matter signature from the rotation curves entirely.

Our analysis yields unphysical dark matter fractions for $\sim 3\%$ of the
sample due to $M_{\mathrm{dyn}}<M_{\mathrm{bar}}$, and an additional 5%
exhibits low $f_{{}_{\mathrm{DM}}}$ (¡0.2) within $R_{opt}$ and $R_{out}$, see
Figure 9 and 20. In order to understand these objects, we carefully examined
their rotation curves, moment maps, and high-resolution images, yet no
apparent anomalies were found. Consequently, our initial guess is that the
stellar masses of these galaxies are likely overestimated. However, it is
important to remark that issues regarding $M_{\mathrm{dyn}}<M_{\mathrm{star}}$
have been previously reported in studies concerning high-$z\ $galaxies,
including works by Wisnioski et al. (2015); Förster Schreiber & Wuyts (2020),
and Sharma et al. (2021a). We propose several conditions under which the
dynamical mass could be smaller than the baryonic mass:

  * •

Incomplete sampling of the velocity distribution: If the velocity measurements
used to estimate the dynamical mass are limited to only a small fraction of
the stars or gas in a galaxy, the resulting estimate of the dynamical mass
could be lower than the actual value, as hinted in Puglisi et al. (2023). In
this analysis, we did our best to extract the the best possible kinematics
information available in our datasets. However, with better data these
estimates will further improve.

  * •

Non-equilibrium conditions: If a galaxy has recently undergone a major
disturbance, such as a merger or collision, the velocity distribution of its
stars and gas may not yet have settled into a stable equilibrium (Riechers et
al., 2014; Lemaux, 2017; Förster Schreiber & Wuyts, 2020, and reference
therein). In this case, the dynamical mass estimate could be lower than the
actual value. In our sample, we visually inspect the galaxy rotation curves
moment maps, and high resolution photometry, and discard all the potential
mergers. However, it is worth noting that our understanding of the merging
history of galaxies at high-$z$ is still limited.

  * •

Non-circular orbits: If the stars and gas in a galaxy are on non-circular
orbits, the velocity measurements used to estimate the dynamical mass will
only provide a lower limit on the actual mass (Binney & Tremaine, 1987;
Cappellari, 2008, references therein). This can only be understood with future
high-resolution spectroscopy. That is, there is room for improvement in our
dynamical masses.

On the contrary, the most massive galaxies in our sample
($\log(M_{\mathrm{star}}\ [\mathrm{M_{\odot}}])>10.7$) exhibit a significant
abundance of dark matter throughout the galactic radius, surpassing the
average dark matter fraction observed in overall sample. Remarkably, the dark
matter halos of these galaxies are heavier than those found in the most
massive systems in the present-day Universe, such as the Milky Way and
Andromeda galaxies (see Figure 8 and Figure 9). These intriguing findings
propose two plausible scenarios:

  1. 1.

These massive galaxies may have followed distinct evolutionary pathways,
resulting in the formation of exceptionally massive dark matter halos. Or,

  2. 2.

The baryonic masses of these objects are underestimated, for instance due to
the presence of massive compact objects. Indeed, if a galaxy contains a
substantial number of compact objects, such as black holes or neutron stars,
their gravitational influence could dominate the motion of stars and gas,
leading to an overestimate of the dynamical mass (Naab & Ostriker, 2017, see
references therein).

Another potential source of uncertainty could arise from variations in stellar
mass estimates due to the use of different SED fitting codes. However, as
demonstrated in Figure 3, the offset from main-sequence (an indirect
measurement of star-formation rate) for all samples, at a fixed stellar mass,
falls within a range of 0.2 dex. Moreover, Figure 21 shows that the stellar
mass within $R_{out}$ across different samples varies by no more than 0.03
dex. This strongly suggests that our results are not dominated by systematic
differences between the SED codes.

Finally, it is important to highlight that galaxies with low baryon surface
densities ($\log(\Sigma_{\text{bar}}(<R_{\text{out}}))<8.0$), often
categorized as low surface-brightness galaxies, are noticeably absent from our
sample at $z>1$, as depicted in Figure 12 (bottom left panel). This absence is
potentially attributable to Tolman Dimming, a relativistic phenomenon
described by general relativity in which the observed surface brightness of a
celestial object diminishes as $(1+z)^{2}$ (Tolman, 1930; Pahre et al., 1996;
Sandage, 2010). Consequently, low surface-brightness galaxies could go missing
in high-redshift observations. It should also be noted that these galaxies
present challenges for observation even in the local Universe. To acquire
resolved rotation curves for such systems, observatories with substantial
aperture sizes, such as the forthcoming Extremely Large Telescope, would be
required.

With these caveats in mind, it is noteworthy that the majority (approximately
75%) of our sample consists of dark matter-dominated galaxies, consistent with
the findings of a previous study conducted by Sharma et al. (2021b), and not
very different from local galaxies. Consequently, the reliability of our
sample, techniques, and measurements instills confidence and provides a firm
foundation for further discussion regarding the evolution of dark matter in
galaxies across cosmic time.

##  6 Discussion

In this section, we discuss the main findings of this work.

Dark Matter Fraction: In local star-forming galaxies, where dark matter is
believed to constitute the majority of the mass in most galaxies, dark matter
fraction estimates typically range from 70% to 90% of the total mass. The dark
matter fraction in these galaxies is relatively higher in the outskirts
compared to the inner regions, indicating that the inner regions are dominated
by baryonic processes while the dynamics of the outskirts is governed by dark
matter (Rubin et al., 1980; Persic et al., 1996; Martinsson et al., 2013;
Courteau & Dutton, 2015; McGaugh, 2016). A similar trend is observed in
high-$z\ $ galaxies, as shown in Figure 9 & 10. Based on the new datasets
presented in this study, we show that the galaxies at high redshifts (0.5 ¡ z
¡ 2.5) are predominantly influenced by dark matter, with fractions ranging
from 50% to 90%. On average, the dark matter fraction appears to decrease as a
function of redshift between $0\leq z\leq 1.8$ across all galactic scales, but
it does not fall below 50%. Measurements at higher redshifts ($z>1.8$) are
currently not reliable due to the limited/small sample size.

While we identify dark matter-dominated systems at high-$z\ $, we observe a
significant scatter in relations, such as $f_{{}_{DM}}-M_{star}$,
$f_{{}_{DM}}-z$, $f_{{}_{DM}}-V_{Rout}$, and $\Sigma_{bar}-f_{{}_{DM}}$, see
Figure 9, 10, & 11, respectively. A similar scatter is also reported in Sharma
et al. (2021b). The scatters in these relations suggest that these galaxies
might still be undergoing the process of building or acclimating their
distribution of baryons and dark matter. In other words, they are at different
stages of ‘galaxy assembly’.121212Galaxy assembly refers to the process by
which galaxies form and evolve over time, involving the accretion of gas and
dark matter, enhanc star formation, and the merging of smaller galaxies into
larger structures. The scatter in the dark matter fraction itself may stem
from various other factors, including differences in formation history of
galaxies, the irregular distribution of baryons that can affect the
distribution of dark matter, and the environment in which these galaxies
reside (Dutton et al., 2016; Behroozi et al., 2019). Moreover, it could also
be due to the diversity in the dark matter halo properties (e.g.,
concentration and expansion), which are closely coupled to the properties of
the baryonic matter (e.g., halo spin parameter and baryon fraction). However,
note that if any systematic uncertainties in estimating baryonic masses have
been neglected, this would naturally increase the scatter, never decrease it,
i.e., reported scatter is a clear upper limit.

Dark Matter Halo Assembly: According to the current cosmological model, Lambda
Cold Dark Matter ($\Lambda$CDM), dark matter plays a fundamental role in the
assembly of cosmic structures in the Universe. This assembly process follows a
hierarchical pattern, wherein smaller structures form first and subsequently
merge to form larger ones, and gradually increase in mass and size (Peebles,
1993). This hierarchical merging process is primarily driven by gravity and
influenced by the distribution of matter in the Universe, encompassing both
dark and baryonic matter. Since dark matter is invisible and cannot be
directly observed, studying the distribution and properties of baryons, such
as surface density and motion governed by the total gravitational potential,
provides a way to gain insights into the distribution of dark matter and its
assembly history.

In accordance with the aforementioned concepts, we investigated
$\Sigma_{bar}-f_{{}_{DM}}$ relation, which was previously examined by Genzel
et al. (2020) within the effective radius. Expanding upon their work, we
extended this relation to encompass the outer radius. Notably, our data well
fit the relationship described by Equation 4, with an intrinsic scatter of
$0.13$ dex. Our best-fit result is in agreement with the relation presented by
Genzel et al. (2020) within a $1\sigma$ uncertainty, as shown in the left
panel of Figure 11. Intriguingly, we observed that the
$\Sigma_{bar}-f_{{}_{DM}}$ relation maintains a consistent slope also when
examined at the optical radius, and when examined at the outer radius in both
low and high mass systems, as illustrated in the upper panel of Figure 12.
Moreover, this relationship remains unchanged across different redshift
ranges, as demonstrated in bottom right panel of Figure 12. Thus, the
$\Sigma_{bar}-f_{{}_{DM}}$ relation exhibits a universal nature, implying that
the influence of baryonic matter on dark matter appears to be similar across
various stellar mass ranges, and redshift intervals. Higher quality data will
allow to pinpoint the exact scatter of this relation: taking the current
scatter at face value would indeed indicate that the relation is not
fundamental (Milgrom, 1983; Famaey & McGaugh, 2012) but that galaxies are
likely at distinct evolutionary stages, for instance in the process of
establishing their respective disks.

On the other hand, when analyzing the relationship between
$V_{c}-f_{{}_{DM}}$, as depicted in the right panel of Figure 11, the
correlation was a bit tighter, with an intrinsic scatter of $0.11$ dex.
However, this relation is not at all universal, as further exploration of this
relationship for low and high stellar mass galaxies, as illustrated in Figure
13, revealed the emergence of distinct sequences as a function of both mass
and redshift. Galaxies with high stellar mass exhibit a lower fraction of dark
matter at fixed rotational velocity. In the lower panel of Figure 13, we
divide the $f_{{}_{DM}}-V_{c}$ relation of low and high stellar masses in low
and high-z ($z\leq 1$ and $z>1$, respectively). Notably, a distinct offset was
observed between $z\leq 1$ and $z>1$ objects for the low-mass bin
($\log(M_{star}\ [M_{\odot}])\leq 10^{10}$). Specifically, low-mass galaxies
at high redshift $z>1$ exhibit a lower fraction of dark matter at fixed
rotational velocity. This is intriguing, suggesting that low-mass mass
galaxies most likely undergo higher degree of evolution in terms of their
respective dark matter and baryon distribution. On the other hand, high-mass
systems already seems to be settled at higer-$z$, as evidenced in the lower
right panel of Figure 13.

##  7 Summary & Conclusions

This study aims at investigating the fraction of dark matter across different
galactic scales and cosmic time. To achieve this objective, we utilized a
substantial sample consisting of 263 main-sequence star-forming disk-like
galaxies within the redshift range $0.6\leq z\leq 2.2$. This sample
encompasses 73 galaxies from the KMOS3D survey, 21 from the KGES survey, and
169 from the KROSS survey. We performed the 3D-forward modeling of datacubes
using the 3DBarolo code, as described in Section 3. Figure 5 provides an
illustrative example of the kinematic modeling results obtained using
3DBarolo, displaying the velocity distribution (moment-1): data, model, and
residuals, followed by the major axis PV-diagram and the velocity dispersion
curve. We applied pressure support corrections to the rotation curves inferred
from 3DBarolo, employing the same methodology established in Sharma et al.
(2021a). These rotation curves are corrected for beamsmearing, inclination,
and pressure gradients, and are referred to as intrinsic rotation curves.

To estimate the dark matter fraction, we subtract the baryon contribution from
the intrinsic rotation curves, assuming that both star and gas are distributed
in an exponential thin disk (Freeman, 1970). The residual contribution is
assumed to originate from the dark matter, required to explain the observed
kinematics. For a detailed explanation of this methodology, we refer the
reader to Section B, C, and Sharma et al. (2021b). We estimated the dark
matter fraction within $R_{e}\ \mathrm{till}\ R_{out}$, see Section 4.1 & 4.2,
and Figure 8 & 9. Additionally, we explored the variation of the dark matter
fraction across different redshift ranges, presented in Section 4.3 and
illustrated in Figure 10. Furthermore, we investigated the interplay between
baryonic matter and dark matter by examining the scaling relations of dark
matter with baryon surface density and circular velocity at the outer radius
of galaxies in Section 4.4, as shown in Figure 11, 12, and 13. The results
obtained from these analyses offer valuable insights into the role of dark
matter in shaping galaxies. Our key findings are:

  * •

Similar to local disk galaxies, high-$z\ $galaxies also exhibit dark matter
dominated outer halos, with dark matter fraction ($f_{{}_{\rm DM}}$) ranging
between 50% and 90%.

  * •

We provide novel observational evidence that the dark matter fraction in
high-$z\ $disk-like galaxies gradually increases from their inner regions to
outskirts, similar to their local counterparts.

  * •

The dark matter fraction demonstrates slightly decreasing trend with
increasing redshifts, as shown in Figure 10. However, it is noteworthy that
even at higher redshifts, the median value of $f_{{}_{\rm DM}}$ remains above
50% across all the galactic scales, very similar to local disk galaxies.

  * •

The $V_{c}-f_{{}_{\rm DM}}$ relation reveals distinct sequences with stellar
mass and redshift intervals. Specifically, low stellar mass ($M_{star}\leq
10^{10}\ {\rm M_{\odot}}$) galaxies exhibit a higher degree of dynamical
evolution, with clearly lower dark matter fractions at a given rotational
velocity at higher redshift ($z>1$). This can potentially be attributed to the
hierarchical assembly of dark matter halos, governing the formation and
evolution of galaxies.

  * •

The $\Sigma_{bar}-f_{{}_{\rm DM}}$ relation demonstrates a consistent trend,
slope, and scatter across different stellar mass ranges and redshift
intervals, indicating a universal nature of the influence of baryonic matter
on the distribution of dark matter. The relation also holds at different radii
within galaxies. This implies that the interplay between baryons and dark
matter in galaxies follows a consistent physics, regardless of their mass,
size, and cosmic time.

Overall, these results represent an important step forward in our
understanding of galaxies and dark matter halo evolution. However, to
establish tighter constraints on these findings, perhaps, we need amendment in
our current kinematic modelling techniques, which heavily rely on the
assumptions of axisymmetry and dynamical equilibrium, that might not be the
case at high-$z$. It is probable that low-surface-brightness galaxies are
under-represented at high redshifts, most-likely falling into the low signal-
to-noise domain and consequently being discarded; hence, making high-quality
data indispensable for resolving the present conundrums. In conclusion,
further exploration and investigation are necessary to unravel the
complexities of dark matter’s influence on baryons, and vice versa, in the
evolution of galaxies.

### Acknowledgments

G.S., acknowledges SARAO postdoctoral fellowship (UId No.: 97882), and thanks
Ambica Govind for providing HST images of full sample. G.S. also thanks Mihael
Petac for various fruitful discussions. J.F. and B.F. acknowledge funding from
the European Research Council (ERC) under the European Unions Horizon 2020
research and innovation program (grant agreement No. 834148) and from the
Agence Nationale de la Recherche (ANR projects ANR18-CE31-0006 and
ANR-19-CE31-0017). GvdV acknowledges funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme under grant agreement No 724857 (Consolidator GrantArcheoDyn). F.R.,
acknowledges support from the Knut and Alice Wallenberg Foundation, and from
the University of Strasbourg Institute for Advanced Study (USIAS), within the
French national programme Investment for the Future (Excellence Initiative)
IdEx-Unistra. M.M., thanks financial support of the Flemish Fund for
Scientific Research (FWO-Vlaanderen), research project G030319N.

## References

  * Angulo & Hahn (2022) Angulo, R. E. & Hahn, O. 2022, Living Reviews in Computational Astrophysics, 8, 1
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143
  * Binney & Tremaine (1987) Binney, J. & Tremaine, S. 1987, Galactic dynamics
  * Boselli et al. (2002) Boselli, A., Lequeux, J., & Gavazzi, G. 2002, A&A, 384, 33
  * Bosma (1981) Bosma, A. 1981, The Astronomical Journal, 86, 1791
  * Bouché et al. (2022) Bouché, N. F., Bera, S., Krajnović, D., et al. 2022, A&A, 658, A76
  * Brammer et al. (2012) Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012, ApJS, 200, 13
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Bullock & Boylan-Kolchin (2017) Bullock, J. S. & Boylan-Kolchin, M. 2017, ARAA, 55, 343
  * Burkert (1995) Burkert, A. 1995, apjl, 447, L25
  * Burkert et al. (2010) Burkert, A., Genzel, R., Bouché, N., et al. 2010, The Astrophysical Journal, 725, 2324
  * Burkert et al. (2016) Burkert, A., Schreiber, N. F., Genzel, R., et al. 2016, The Astrophysical Journal, 826, 214
  * Calzetti et al. (1994) Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582
  * Cappellari (2008) Cappellari, M. 2008, MNRAS, 390, 71
  * Cappellari et al. (2013) Cappellari, M., Scott, N., Alatalo, K., et al. 2013, MNRAS, 432, 1709
  * Casey et al. (2023) Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31
  * Chabrier (2003) Chabrier, G. 2003, Publications of the Astronomical Society of the Pacific, 115, 763
  * Charlot & Fall (2000) Charlot, S. & Fall, S. M. 2000, ApJ, 539, 718
  * Chowdhury et al. (2022) Chowdhury, A., Kanekar, N., & Chengalur, J. N. 2022, ApJL, 941, L6
  * Cormier et al. (2016) Cormier, D., Bigiel, F., Wang, J., et al. 2016, Monthly Notices of the Royal Astronomical Society, 463, 1724
  * Courteau & Dutton (2015) Courteau, S. & Dutton, A. A. 2015, ApJL, 801, L20
  * da Cunha et al. (2008) da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595
  * Daddi et al. (2010) Daddi, E., Elbaz, D., Walter, F., et al. 2010, ApJL, 714, L118
  * Dekel et al. (2021) Dekel, A., Freundlich, J., Jiang, F., et al. 2021, MNRAS, 508, 999
  * Di Teodoro et al. (2016) Di Teodoro, E., Fraternali, F., & Miller, S. 2016, Astronomy & Astrophysics, 594, A77
  * Dutton et al. (2016) Dutton, A. A., Macciò, A. V., Dekel, A., et al. 2016, MNRAS, 461, 2658
  * Famaey & McGaugh (2012) Famaey, B. & McGaugh, S. S. 2012, Living Reviews in Relativity, 15, 10
  * Förster Schreiber et al. (2006) Förster Schreiber, N. M., Genzel, R., Lehnert, M. D., et al. 2006, ApJ, 645, 1062
  * Freeman (1970) Freeman, K. 1970, The Astrophysical Journal, 160, 811, provided by the SAO/NASA Astrophysics Data System
  * Freundlich et al. (2019) Freundlich, J., Combes, F., Tacconi, L. J., et al. 2019, A&A, 622, A105
  * Freundlich et al. (2020) Freundlich, J., Dekel, A., Jiang, F., et al. 2020, MNRAS, 491, 4523
  * Fu et al. (2010) Fu, J., Guo, Q., Kauffmann, G., & Krumholz, M. R. 2010, Monthly Notices of the Royal Astronomical Society, 409, 515
  * Förster Schreiber & Wuyts (2020) Förster Schreiber, N. M. & Wuyts, S. 2020, Annual Review of Astronomy and Astrophysics, 58, null
  * Genzel et al. (2020) Genzel, R., Price, S. H., Übler, H., et al. 2020, ApJ, 902, 98
  * Genzel et al. (2017) Genzel, R., Schreiber, N. F., Ubler, H., et al. 2017, Nature, 543, 397
  * Genzel et al. (2015) Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20
  * Giacconi et al. (2001) Giacconi, R., Rosati, P., Tozzi, P., et al. 2001, The Astrophysical Journal, 551, 624
  * Gillman et al. (2020) Gillman, S., Tiley, A. L., Swinbank, A. M., et al. 2020, MNRAS, 492, 1492
  * Harris et al. (2020) Harris, W. E., Remus, R.-S., Harris, G. L. H., & Babyk, I. V. 2020, ApJ, 905, 28
  * Harrison et al. (2017) Harrison, C., Johnson, H., Swinbank, A., et al. 2017, Monthly Notices of the Royal Astronomical Society, 467, 1965
  * Häußler et al. (2013) Häußler, B., Bamford, S. P., Vika, M., et al. 2013, MNRAS, 430, 330
  * Johnson et al. (2018) Johnson, H. L., Harrison, C. M., Swinbank, A. M., et al. 2018, Monthly Notices of the Royal Astronomical Society, 474, 5076
  * Kassin et al. (2006) Kassin, S. A., de Jong, R. S., & Weiner, B. J. 2006, ApJ, 643, 804
  * Keres et al. (2003) Keres, D., Yun, M. S., & Young, J. S. 2003, ApJ, 582, 659
  * Kriek et al. (2009) Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221
  * Lagos et al. (2011) Lagos, C. D. P., Baugh, C. M., Lacey, C. G., et al. 2011, MNRAS, 418, 1649
  * Lang et al. (2017) Lang, P., Schreiber, N. M. F., Genzel, R., et al. 2017, The Astrophysical Journal, 840, 92
  * Lawrence et al. (2007) Lawrence, A., Warren, S., Almaini, O., et al. 2007, Monthly Notices of the Royal Astronomical Society, 379, 1599
  * Lehmer et al. (2005) Lehmer, B. D., Brandt, W., Alexander, D., et al. 2005, The Astrophysical Journal Supplement Series, 161, 21
  * Lelli et al. (2023) Lelli, F., Zhang, Z.-Y., Bisbas, T. G., et al. 2023, arXiv e-prints, arXiv:2302.00030
  * Lemaux (2017) Lemaux, B. C. 2017, in Early stages of Galaxy Cluster Formation, 43
  * Leroy et al. (2008) Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782
  * Li et al. (2023) Li, Z., Dekel, A., Mandelker, N., Freundlich, J., & François, T. L. 2023, MNRAS, 518, 5356
  * Martin et al. (2010) Martin, A. M., Papastergis, E., Giovanelli, R., et al. 2010, apj, 723, 1359
  * Martinsson et al. (2013) Martinsson, T. P. K., Verheijen, M. A. W., Westfall, K. B., et al. 2013, A&A, 557, A131
  * Martorano et al. (2023) Martorano, M., van der Wel, A., Bell, E. F., et al. 2023, arXiv e-prints, arXiv:2308.11392
  * McGaugh (2010) McGaugh, S. 2010, in American Institute of Physics Conference Series, Vol. 1240, Hunting for the Dark: the Hidden Side of Galaxy Formation, ed. V. P. Debattista & C. C. Popescu, 13–16
  * McGaugh (2016) McGaugh, S. S. 2016, ApJ, 816, 42
  * Milgrom (1983) Milgrom, M. 1983, ApJ, 270, 365
  * Momcheva et al. (2016) Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27
  * Naab & Ostriker (2017) Naab, T. & Ostriker, J. P. 2017, ARAA, 55, 59
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, apj, 462, 563
  * Nestor Shachar et al. (2022) Nestor Shachar, A., Price, S. H., Förster Schreiber, N. M., et al. 2022, arXiv e-prints, arXiv:2209.12199
  * Noeske et al. (2007) Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJL, 660, L43
  * Pahre et al. (1996) Pahre, M. A., Djorgovski, S. G., & de Carvalho, R. R. 1996, ApJL, 456, L79
  * Peebles (1993) Peebles, P. J. E. 1993, Principles of Physical Cosmology
  * Peng et al. (2010) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097
  * Persic & Salucci (1990) Persic, M. & Salucci, P. 1990, mnras, 247, 349
  * Persic et al. (1996) Persic, M., Salucci, P., & Stel, F. 1996, Monthly Notices of the Royal Astronomical Society, 281, 27
  * Petač (2020) Petač, M. 2020, Phys. Rev. D, 102, 123028
  * Pontzen & Governato (2012) Pontzen, A. & Governato, F. 2012, MNRAS, 421, 3464
  * Pontzen & Governato (2014) Pontzen, A. & Governato, F. 2014, Nature, 506, 171
  * Price et al. (2016) Price, S. H., Kriek, M., Shapley, A. E., et al. 2016, ApJ, 819, 80
  * Puglisi et al. (2023) Puglisi, A., Dudzevičiūtė, U., Swinbank, M., et al. 2023, arXiv e-prints, arXiv:2305.04382
  * Renaud et al. (2021) Renaud, F., Romeo, A. B., & Agertz, O. 2021, MNRAS, 508, 352
  * Riechers et al. (2014) Riechers, D. A., Carilli, C. L., Capak, P. L., et al. 2014, ApJ, 796, 84
  * Rodighiero et al. (2011) Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJL, 739, L40
  * Rubin et al. (1980) Rubin, V. C., Ford Jr, W. K., & Thonnard, N. 1980, The Astrophysical Journal, 238, 471
  * Salucci (2019) Salucci, P. 2019, A&AR, 27, 2
  * Sandage (2010) Sandage, A. 2010, AJ, 139, 728
  * Scoville et al. (2007) Scoville, N., Abraham, R. G., Aussel, H., et al. 2007, The Astrophysical Journal Supplement Series, 172, 38
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Sharma et al. (2021a) Sharma, G., Salucci, P., Harrison, C. M., van de Ven, G., & Lapi, A. 2021a, MNRAS, 503, 1753
  * Sharma et al. (2021b) Sharma, G., Salucci, P., & van de Ven, G. 2021b, A&A, 653, A20
  * Sharma et al. (2022) Sharma, G., Salucci, P., & van de Ven, G. 2022, A&A, 659, A40
  * Sharples (2014) Sharples, R. 2014, Proceedings of the International Astronomical Union, 10, 11–16
  * Skelton et al. (2014) Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24
  * Sobral et al. (2015) Sobral, D., Matthee, J., Best, P. N., et al. 2015, Monthly Notices of the Royal Astronomical Society, 451, 2303
  * Sofue & Rubin (2001) Sofue, Y. & Rubin, V. 2001, ARAA, 39, 137
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, The Astrophysical Journal Supplement Series, 214, 15
  * Steidel et al. (1998) Steidel, C. C., Adelberger, K. L., Dickinson, M., et al. 1998, The Astrophysical Journal, 492, 428
  * Stott et al. (2016) Stott, J. P., Swinbank, A., Johnson, H. L., et al. 2016, Monthly Notices of the Royal Astronomical Society, 457, 1888
  * Tacconi et al. (2018) Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179
  * Tamm et al. (2012) Tamm, A., Tempel, E., Tenjes, P., Tihhonova, O., & Tuvikene, T. 2012, A&A, 546, A4
  * Teodoro & Fraternali (2015) Teodoro, E. D. & Fraternali, F. 2015, Monthly Notices of the Royal Astronomical Society, 451, 3021
  * Tiley et al. (2019a) Tiley, A., Bureau, M., Cortese, L., et al. 2019a, Monthly Notices of the Royal Astronomical Society, 482, 2166
  * Tiley et al. (2021) Tiley, A. L., Gillman, S., Cortese, L., et al. 2021, MNRAS, 506, 323
  * Tiley et al. (2019b) Tiley, A. L., Swinbank, A., Harrison, C., et al. 2019b, Monthly Notices of the Royal Astronomical Society, 485, 934
  * Tolman (1930) Tolman, R. C. 1930, Proceedings of the National Academy of Science, 16, 511
  * Tully & Fisher (1977) Tully, R. B. & Fisher, J. R. 1977, A&A, 54, 661
  * Übler et al. (2019) Übler, H., Genzel, R., Wisnioski, E., et al. 2019, ApJ, 880, 48
  * Valentino et al. (2023) Valentino, F., Brammer, G., Gould, K. M. L., et al. 2023, ApJ, 947, 20
  * van Albada et al. (1985) van Albada, T. S., Bahcall, J. N., Begeman, K., & Sancisi, R. 1985, ApJ, 295, 305
  * van Dokkum & Brammer (2010) van Dokkum, P. G. & Brammer, G. 2010, ApJL, 718, L73
  * Weijmans et al. (2014) Weijmans, A.-M., de Zeeuw, P. T., Emsellem, E., et al. 2014, MNRAS, 444, 3340
  * Whitaker et al. (2012) Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJL, 754, L29
  * Wisnioski et al. (2019) Wisnioski, E., Förster Schreiber, N. M., Fossati, M., et al. 2019, ApJ, 886, 124
  * Wisnioski et al. (2015) Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209
  * Wuyts et al. (2011) Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011, ApJ, 742, 96
  * Wuyts et al. (2016) Wuyts, S., Förster Schreiber, N. M., Wisnioski, E., et al. 2016, ApJ, 831, 149
  * Zwaan et al. (2005) Zwaan, M. A., Meyer, M. J., Staveley-Smith, L., & Webster, R. L. 2005, Monthly Notices of the Royal Astronomical Society: Letters, 359, L30
  * Übler et al. (2018) Übler, H., Genzel, R., Tacconi, L. J., et al. 2018, The Astrophysical Journal, 854, L24

##  Appendix A Constraining Gas Geometrical Parameters with 3DBarolo

Firstly, observations of high-$z\ $galaxies are often poorly resolved and
feature low signal-to-noise ratios (SN), making kinematic modeling a
challenging task. Therefore, in order to obtain the robust results, we
investigate their SN and H$\alpha$ images, as shown in Figures 14. Based on
this quantitative and qualitative selection, we divide the sample in three
categories: Q1: SN$>3$ and sharp H$\alpha$ image; Q2: SN$\geq 3$ and
moderately visible source; Q3: either SN$<3$ or no appearance of source. We
discard all the Q3 galaxies.

Secondly, as we go higher in redshift the information on morphological gas
central x-y position ($x_{c},y_{c}$) and position angle (PA) starts
differentiating from their photometric measurements (as reported in Wisnioski
et al. 2015; Harrison et al. 2017; Sharma et al. 2021a). On the other hand,
from our previous work with 3DBarolo (Sharma et al. 2021a), we learned that
3DFIT TASK is not feasible in constraining the multiple parameters. Therefore,
in this work we estimate gas geometrical parameters using an optimization
function that runs atop of 3DBarolo run. This optimization function takes
following inputs: (1) loss function given in Eq. 3, (2) Initial guess and
bounds on parameter to be fitted, (3) minimization method, and (4) ‘par file’
to execute 3DBarolo run. In particular, we optimize PA, $x_{c}$, and $y_{c}$ .
Initial guess on these parameters comes from the photometric information or
from H$\alpha$ datacubes (discussed in Sec. 2). The bounds are as following:
$[x_{c}=x_{phot}\pm 4,y_{c}=y_{phot}\pm 4]$ and PA[0, 360]. To be accurate, we
also optimize systemic velocity ($V_{sys}$), the initial gas on it is computed
from ‘Doppler shift’ of H$\alpha$ line, and bounds are [-200, 200]. However,
choice of fitting $V_{sys}$ is tricky, for example, when we observe high
residuals in moment maps, we fix it to a median value given in 3DBarolo output
files. The optimize-minimize function works in three steps:

  1. 1.

Run 3D Barolo (free: $V_{rot}$, $V_{disp}$, and in some cases $V_{sys}$),
which gives moment maps, PV- diagrams, and estimates on $V_{rot}$ and
$V_{disp}$.

  2. 2.

Compute the loss function (Eq. 3). The data and model inputs in loss function
depend on the choice of parameter to be fitted, which we have discussed below.

  3. 3.

Save the loss, and optimization log. The optimization log gives us the best
estimate of parameters and success rate (True or False).

Figure 14: An example of quality assessment of KMOS3D and KGES sample. Row 1-2
shows the Q1 and Q2 galaxies. Row 3-4 are Q3 galaxies. The color codes are
given in the legend of the plots. Figure 15: This figure compares the circular
velocity of KROSS galaxies as computed in this study (y-axis) with that of
Sharma et al. (2021b, x-axis). The Quality-1 and Quality-2 galaxies, as
classified in Sharma et al. (2021a), are represented by yellow and purple
circles, respectively, and the black line shows one-to-one relationship
between two studies.

The optimization procedure follows a logical flow consists of four stages.
Firstly, the $x_{c}$ and $y_{c}$ parameters are constrained using moment-0
maps, which provide information on the gas distribution. In the second stage,
the previously optimized $x_{c}$ and $y_{c}$ positions are fixed to optimize
the position angle (PA) parameter using minor axis PV-diagrams. As noted in
(Sharma et al. 2021a, appendix-A), the symmetry of the minor axis PV-diagram
around the axis indicates the correctness of the PA parameter; an asymmetric
PV-diagram implies an incorrect PA. In the third stage, the previously
constrained $x_{c}$, $y_{c}$, and PA parameters are fixed to estimate the
systemic velocity ($V_{sys}$) parameter. The loss function for $V_{sys}$ uses
major axis PV-diagrams since $V_{sys}$ directly affects the rotation velocity,
which is determined by the major axis PV-diagram. In the final stage, all
optimized parameters are used to compute the final loss using moment-1 maps
since our goal is to constrain the rotation velocity of the system. However,
for cross-check, the final loss is also computed using moment-2 maps, although
the difference is negligible ($<0.001$) and does not appear in the residual
maps. Lastly, the output corresponding to the minimum loss is saved.

We noted that this optimization procedure does not work for random guesses on
the parameter values and requires reasonably close initial estimates, such as
those provided by photometric estimates. Additionally, the procedure is
computationally intensive, with each galaxy taking approximately 20-30 minutes
to optimize when running on 8-cores. Nevertheless, this procedure is
worthwhile for studying the nature of dark matter.

We evaluated the effectiveness of 3DBarolo+optimizer on the KROSS data, which
we previously analyzed and discussed in Sharma et al. (2021a). As expected,
the optimization process failed for galaxies with low SN ($<5$), specifically
those with Quality2. Consequently, we were compelled to exclude 56 out of 225
galaxies. In Figure 15, we compare the circular velocities (computed at
$R_{out}$) of KROSS sample, obtained with old and new kinematic modeling
techniques. It is noteworthy that the previous and current circular velocity
measurements are in reasonable agreement, with an intrinsic scatter of 0.03
dex. Both Quality-1 and Quality-2 galaxies in the GS21 sample exhibit
variations in their circular velocities. Nonetheless, we used the new
measurements because assessment of data-to-model fits in the position
velocity-diagram (PV-diagram) are better in the new analysis, as shown in
Appendix Figure LABEL:fig:KROSS-BB (Fig. is available upon request or see A&A
version of Article). Moreover, moment-map residuals stays close to zero.
Therefore, we conclude that 3DBarolo+optimizer works well on high-$z$ moderate
signal-to-noise data, and hence can be used to constrain extra parameters.

Figure 16: This figure shows the rotation to dispersion ratio ($v/\sigma$) of
galaxies as a function of circular velocity computed at $R_{out}$. The left
and right panels show $v/\sigma$ before and after pressure gradient
corrections, respectively. The red horizontal dashed line indicates
$v/\sigma=1.1$. This figure demonstrates that, after pressure support
corrections, all the galaxies are rotation-dominated. Therefore, we do not
exclude any galaxies based on their intrinsic $v/\sigma$ values.

Figure 17: A few examples of mass modeled rotation curves in case of maximum
stellar disk scenario (upper panel) and maximum stellar+gas disk (lower
panel). The best best fit parameters are shown in right column, and color
codes are given in the legend of the plot. This example is favouring the
maximum stellar and baryonic disk scenario, upper and lower panels,
respectively.

Figure 18: A few examples of mass modeled rotation curves in case of maximum
stellar disk scenario (upper panel) and maximum stellar+gas disk (lower
panel). The best best fit parameters are shown in right column, and color
codes are given in the legend of the plot. This example demonstrates that (1)
the model is incapable of fitting the observations using only baryons, and (2)
the best-fit values for stellar or gas mass are an order of magnitude higher
than their fiducial values, which is unrealistic given the uncertainties in
the observations. Therefore, we suggest that the dark matter halo component is
required in high-$z\ $galaxies.

##  Appendix B Mass modeling of rotation curves

Initially, we assume the absence of dark matter halos around high-$z$ galaxies
and dynamically model the baryonic component to fully fit the observed
rotation curves. In this case, observed velocity is defined as:

$V^{2}_{c}(R)=V^{2}_{\mathrm{stars}}(R)+V^{2}_{\mathrm{gas}}(R)$ (6)

where, $V_{\mathrm{stars}}$ represents the contribution of stars in both the
bulge and disk, while $V_{\mathrm{gas}}$ encompasses the combined effect of
molecular and atomic gas components. It is important to note that in the
maximum stellar-disk scenario, the second term on the right-hand side of the
equation becomes null.

We perform the dynamical modeling using a Bayesian inference technique
established in Sharma et al. (2022), which employs flexible models. For
Bayesian inference, we employ Markov chain Monte Carlo (MCMC) analysis,
assuming a $\chi^{2}$ test statistic on observed kinematics and modeled
kinematics. In our modeling approach, we assume that the distribution of stars
and gas follows an exponential mass profile (Freeman 1970), while the bulge is
treated as a point mass. It is worth noting that the bulge remains unresolved
in all datasets. Moreover, we have only fitted the outer rotation curves;
therefore, we have fixed the bulge mass at $10^{9}\ M_{\odot}$ as a first-
order approximation. To estimate the molecular and atomic gas masses, we use
observed scaling relations, as described in Section 2.4, where we fix star-
formation rate at its fiducial value. This approach allows us to obtain
realistic gas masses for a given stellar mass and star-formation rate. In the
modeling we keep flat priors on stellar disk mass ($8.0\leq log(M_{star}\
\mathrm{M_{\odot}})\leq 12.0$) and stellar disk radius keeps the Gaussion
priors with 25% uncertainty around true value. A comprehensive discussion of
this mass modeling technique can be found in Sharma et al. (2022, sec: 3.1).

In Figures 17 and 18, we provide a few examples of mass-modeled rotation
curves (rotation curves) to illustrate different scenarios. The upper panel of
Figure 17 depicts the maximum stellar disk scenario, while the lower panel
represents the maximum stars+gas disk scenario. The corresponding best-fit
parameters are displayed in the right column. Figures 17 demonstrates the
robustness of our mass modeling technique. However, Figure 18 reveals the
limitations of using the same models, as we are unable to precisely constrain
the stellar mass, whether considering the maximum stellar disk or the
stars+gas disk scenario. In some cases depicted in the same figure (lower
panel), the mass modeling allows for extremely high gas masses while
suppressing the stellar masses to values lower than their observed photometric
masses. Furthermore, we observe that the majority of the sample exhibits best-
fit stellar masses that are an order of magnitude higher than their fiducial
values (see Figure 7). These discrepancies are both alarming and unrealistic,
indicating the probable necessity of an additional halo component in these
galaxies.

##  Appendix C Halo Model Independent $f_{\mathrm{{}_{DM}}}$

In our previous work Sharma et al. (2021b), we established a halo model
independent framework to study the dark matter fraction of high-z. In this
framework, first we compute the total dynamical mass of galaxy directly from
its rotation curve. The dynamical mass of a galaxy is defined as

$M_{dyn}(<R)=\kappa(\mathrm{{\tiny R}})\frac{V^{2}(R)\ R}{G},$ (7)

where $V(R)$ is the circular velocity computed at radius $R$ and $\kappa$(R)
is the geometric factor accounting for the presence of a stellar disk
alongside the spherical bulge and halo (see, Persic & Salucci 1990). For
$R_{e}$, $R_{opt}$, and $R_{out}$, the values of $\kappa$ are 1.2, 1.05, and
1.0, respectively (see also Sharma et al. 2021b). The bulge mass contribution
within 5 kpc is negligible for local spirals; therefore, we do not model it.
Moreover, the geometric parameter $\kappa$(R) in the dynamical mass
calculation takes into account the distribution of the mass located in the
bulge and in the disk.

Second, we disentangle the contribution of total baryonic mass from dynamical
mass. 131313Total baryonic mass is $M_{bar}=M_{*}+\ M_{H_{2}}+1.33M_{HI}$,
where factor 1.33 accounts the Helium content.  Briefly, we know that our
samples are good representative of main-sequence of star-forming galaxies as
shown in Figure 3, which allow us to compute their molecular and atomic masses
using scaling relations (Tacconi et al. 2018; Chowdhury et al. 2022; Lagos et
al. 2011, respectively). Our goal is to calculate the dark matter fraction of
our galaxies within $R_{e},\ R_{opt}$, and $R_{out}$, and therefore we must
first determine the baryonic masses within these radii. In (Sharma et al.
2021a), we found that rotation curves of $z\sim 1$ galaxies are similar to
local star-forming disk galaxies. This suggests that the radial distribution
of stellar and molecular gas masses within these galaxies can be well
approximated by the Freeman disk (Freeman 1970):

$\begin{array}[]{r@{}l}M_{*}(<R)=M^{tot}_{*}\
\Big{[}1-\Big{(}1+\frac{R}{R_{*}}\Big{)}\ \exp(\frac{-R}{R_{*}})\Big{]}\\\ \\\
M_{H2}(<R)=M^{tot}_{H2}\ \Big{[}1-\Big{(}1+\frac{R}{R_{H2}}\Big{)}\
\exp(\frac{-R}{R_{H2}})\Big{]}\\\ \\\ M_{HI}(<R)=M^{tot}_{HI}\
\Big{[}1-\Big{(}1+\frac{R}{R_{HI}}\Big{)}\
\exp(\frac{-R}{R_{HI}})\Big{]}\end{array}$ (8)

where stars are assumed to be distributed in the stellar disk ($R_{*}=R_{D}$)
known from photometry, discussed in Section [2]. The molecular gas is
generally distributed outward through the stellar disk (up to the length of
the ionized gas $R_{gas}$); therefore, we take $R_{H2}\equiv R_{gas}$. Here,
we estimate the gas scale length $R_{gas}$ by fitting the $H_{\alpha}$ surface
brightness, see Sharma et al. (2021b, see appendix-C). Moreover, studies of
local disk galaxies have shown that the surface brightness of the HI disk is
much more extended than that of the H2 disk (Fu et al. 2010, see their Fig.
5); see also Leroy et al. (2008) and Cormier et al. (2016). Therefore, we
assume $R_{HI}=2\times R_{H2}$, which is a rough estimate, but still
reasonable considering that at high-$z$ no information is available on the
$M_{HI}$ (or $M_{H2}$) surface brightness distribution. Thus, the Equation [8]
allows us to estimate $M_{*}$, $M_{HI}$, and $M_{H2}$ within different radii
($R_{e},\ R_{opt}$, and $R_{out}$).

Given the information on bayonic and dynamical masses, the dark matter
fractions within radius $R$ can be computed as

$f_{DM}(<R)=1-\frac{M_{bar}(<R)}{M_{dyn}(<R)}$ (9)

Therefore, using Equation[9], we computed $f_{DM}$ within $R_{e},\ R_{opt}$,
and $R_{out}$ for all three datasets. We note that, owing to the limited
spatial resolution in our RCs, the measurements of $f_{{}_{\rm DM}}$ within
$R_{e}$ are less accurate than for $R_{opt}$ and $R_{out}$.

The stellar disk in high-$z\ $galaxies is often unresolved; therefore, its is
ambiguous to know whether disk is exponentially thin (Freeman 1970, well
fitted by) or thick (which can be well modeled via S$\acute{e}$rsic profile
Sérsic 1963.) To address this issue, we estimate the stellar mass within a
scale radius (e.g., $R_{out}$) using both the S$\acute{e}$rsic profile (Sérsic
1963) and Freeman disk (Freeman 1970). Let us remark, the S$\acute{e}$rsic
index of galaxies is set to $n=1.35$, which is an average value of our
samples. The results are shown in Figure 21. We observe that the use of the
Freeman disk yields stellar masses that are, on average, 1.020 times higher,
reducing the dark matter fraction by a factor of 1.07. In other words, the
dark matter fraction derived using the Freeman disk will be underestimated by
only 2%. Therefore, stars (and gas) are assumed to be distributed in
accordance with the Freeman disk (Freeman 1970) same as our previous works
(Sharma et al. 2021b, 2022).

##  Appendix D Extra

We had the opportunity to improve the measurement of $R_{e}$ for a subset of
KROSS sample through the latest observations from the James Webb Space
Telescope (JWST). This particular subset of galaxies is located within the
COSMOS field (Skelton et al. 2014), which has recently been observed by the
COSMOS-WEB team (Casey et al. 2023). For these targets JWST/NIRCam sizes in
the filters F115W, F150W, F277W and F444W were recovered via S$\acute{e}$risc
profile fitting using GALFITM (Häußler et al. 2013) with the same methodology
outlined in Martorano et al. (2023) on mosaics publicly available in the Dawn
JWST Archive141414 https://dawn-cph.github.io/dja/index.html(Valentino et al.
2023). In this study, we employ the $R_{e}$ values obtained from the F115W
band, which closely corresponds to the rest-frame $H\alpha$-wavelength
observed at $z\sim 0.85$.

In Figure 19, the left panel illustrates the comparison between the new JWST-
derived $R_{e}$ and the previously measured $R_{e}$ adopted from KROSS parent
catalog. We observe that the JWST $R_{e}$ values are higher than the previous
measurements but remain within a scatter of 0.2 dex. Although the difference
is small, these newly improved $R_{e}$ measurements are crucial for
investigating the dark matter fraction, especially within $R_{e}$. Therefore,
we have estimated the dark matter fraction using the updated $R_{e}$ values,
as shown in the right panel of Figure 19. We observe that galaxies with low
dark matter fractions (¡20%) have been pushed towards higher values. In this
specific sub-sample, none of the galaxies have a dark matter fraction below
$20\%$.

Figure 19: Left panel: Comparison of previously measured $R_{e}$ values of the
KROSS sub-sample with the new JWST $R_{e}$. The solid black line represents
the one-to-one relation, while the dashed black lines indicate a scatter of
0.2 dex. Right panel: Dark matter fraction within $R_{e}$ using JWST
photometry. Objects are color-coded based on their circular velocities within
$R_{e}$.  Figure 20: Dark matter fraction within $R_{opt}$. The KMOS3D, KGES,
and KROSS datasets are represented in red, green, and blue colors,
respectively. For comparison, we include two local massive disk galaxies: the
Milky Way (MW) and Andromeda (M31), represented by the blue and the orange
star, respectively. The solid black dashed line represents the 100% dark
matter regime, while the yellow dashed line represents the baryon-dominated
regime. The gray shaded area indicates the forbidden region, where galaxies
with $M_{dyn}<M_{bar}$ are located.  Figure 21:  Comparison of stellar masses
derived using S$\acute{e}$rsic profile (y-axis) and Freeman disk (x-axis).
Different color code are representing the different scale radii of galaxies
within which stellar masses are estimated (blue: $R_{out}$, orange: $R_{opt}$,
and yellow: $R_{e}$), and black dotted-dashed line represents the one-to-one
relation.

In Figure 20 , we plot the dark matter fraction within $R_{opt}$ ($\sim 2\
R_{e}$), of galaxies as function of stellar mass. Firstly, we observe that the
majority of the sample exhibits dark matter-dominated outer disks, with 16% of
objects showing $0.2\leq f_{{}_{DM}}(<R_{opt})<0.5$, and only 6% of objects
having $f_{{}_{DM}}(<R_{opt})<0.2$ (including objects in the forbidden
region). Secondly, we notice a slightly decreasing trend of $f_{{}_{DM}}$ as a
function of stellar mass, except for a few outliers. Lastly, we compared the
dark matter fraction within $R_{opt}$ and $R_{out}$. Interestingly, we
observed that galaxies generally exhibit slightly lower dark matter fractions
in the inner region ($R_{opt}$) compared to outskirts ($R_{out}$). On average,
dark matter fraction within $R_{opt}$ is about $70\%$, while it is $\sim 80\%$
within $R_{out}$, these finding coincides with the results reported in Sharma
et al. (2021b).

In Figure 22, we present the binned dark matter fraction as a function of
redshift within three different radii, namely $R_{e}$, $R_{\mathrm{opt}}$, and
$R_{\mathrm{out}}$. To facilitate comparisons and provide context, we have
included measurements from local star-forming galaxies (Persic et al. 1996).
The figure illustrates a decreasing trend in the dark matter fraction between
$0\leq z\leq 1.8$ across all radii. Upon comparing our finding with a latest
high-$z\ $ study (Nestor Shachar et al. 2022), especially for dark matter
fraction within $R_{e}$, we noticed major differences. With current data in
hands, we can not fill this huge gap of discrepancy, more high-resolution data
and further investigation is required.

Figure 22: Dark matter fraction across the cosmic time ($0.0\leq z<2.5$). We
show averaged dark matter fraction within $R_{e}$, $R_{opt}$, and $R_{out}$ as
a function of redshift, represented by yellow, orange, and blue lines,
respectively. The bootstrap-errors are indicated by shaded areas color-coded
with the same color as the averaged line. Dark matter fraction of local star-
forming galaxies within $R_{out}$ and $R_{opt}$ are from Persic et al. (1996)
and within $R_{e}$ are taken from disk galaxy survey Courteau & Dutton (2015).
The uncertainties on these measurements are, on average, $\pm 0.25$. The color
code of these estimates are same as high-$z$ galaxies. The dark-yellow dotted-
dashed lines shows the $f_{{}_{DM}}(<R_{e})$ from Nestor Shachar et al.
(2022).

Figure 23: HST Images and RCs of massive galaxies having high dark matter
fraction, i.e., galaxies those are standing odd in $f_{DM}-M_{star}$ relation.