Observational Evidence of Evolving Dark Matter Profiles at z≤ 1

Gauri Sharma,
1,2,3Contact<EMAIL_ADDRESS>Paolo Salucci 1,2,3
Glenn van de Ven4

SISSA International School for Advanced Studies, Via Bonomea 265, I-34136 Trieste, Italy
GSKY, INFN-Sezione di Trieste, via Valerio 2, I-34127 Trieste, Italy
IFPU Institute for Fundamental Physics of the Universe, Via Beirut, 2, 34151 Trieste, Italy
Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17, 1180 Wien, Austria

In the concordance cosmological scenario the cold, collisionless  Dark Matter component dominates the mass budget of galaxies and interacts with baryons only via gravity. However, there is growing evidence, that the former, instead,  responds to the baryonic (feedback) processes by modifying its density distribution. These processes can be captured by comparing the inner dynamics of galaxies across the  cosmic-time.

A pilot study of dynamical mass modelling of high redshift galaxy rotation curves, which is capable of constraining the structure of dark matter halos across cosmic time.

We investigate the dark matter halos of 256 star-forming disc-like galaxies at $z\sim 1$ using the KMOS redshift one spectroscopic survey (KROSS). This sample covers the redshifts $0.6 \leq z \leq 1.04$, effective radii $0.69 \leq R_e [\mathrm{kpc}] \leq 7.76$, and total stellar masses $8.7 \leq log(M_{\mathrm{star}} \ [\mathrm{M_\odot}]) \leq 11.32$. We present a mass modelling approach to study the rotation curves of these galaxies, which allow us to dynamically calculate the physical properties associated with the baryons and the dark matter halo. For the former we assume a Freeman disc, while for the latter we employ the NFW and the Burkert halo profiles, separately. At the end, we compare the results of both cases with state-of-the-art cosmological galaxy simulations (EAGLE, TNG100 and TNG50).

We find that the cored dark matter halo emerged as the  dominant quantity from  a radius 1-3 times the effective radius. Its fraction to the total mass is in good agreement with the outcome of hydrodynamical galaxy simulations. Remarkably, we found  that the dark matter core of $z\sim 1$ star-forming galaxies are smaller and denser than their local counterparts.

We conclude that dark matter halos have gradually expanded over the past 6.5 Gyrs. That is, observations are capable of capturing the dark matter response to the baryonic processes (e.g. feedbacks), and thus giving us the first empirical evidence of gravitational potential fluctuations in the inner region of galaxies, which can be verified with deep surveys and future missions.

§ INTRODUCTION

Dark Matter (DM) is a type of matter that astrophysicists have brought forward to account for the observational effects that appear to be due to the presence of an invisible massive component <cit.>. The astrophysical nature of such invisible matter can be inferred from its gravitational interaction on the luminous matter <cit.> and on the radiation which emerges at different cosmological scales <cit.>. On the other hand, its particle nature, which is the centre of DM phenomenon is not yet characterized.

Recently, DM phenomenon was uniquely framed into the dark energy driven cold dark matter ($\mathrm{\Lambda CDM}$) scenario <cit.>, where non-relativistic DM can be described as a collision-less fluid, whose particle with cross-section of about $10^{-26} \ \mathrm{cm^{2}}$, only interact gravitationally and very weakly with standard model particles.  This scenario, which emerges cosmologically very clearly <cit.>, has been confirmed by large scale observations <cit.> and can be investigated by appropriate simulations <cit.> and semi-analytical modelling <cit.>. As a consequence, a very well defined perturbation power law spectrum and lack of dependencies on the initial conditions, yields to a specific bottom-up formation and evolution of DM halos. However, two aspects in the progressive field of DM opposes the latter: 1) The search of DM particles, performed by various direct and indirect methods <cit.>, one of the major efforts of astro-particle physics, has been so far unsuccessful <cit.>; 2) the main imprint of collision-less particle is a  DM halo with cuspy central DM density ($\rho \propto r^{-1}$), which is not present in the majority of local spirals, dwarf-discs and low surface brightness galaxies <cit.>. Moreover, local galaxies show that the physical properties (e.g., luminosity and disc scale length) of luminous matter are deeply connected with the structural properties of DM (e.g. its core radius and core density), see <cit.> and <cit.>. This is very difficult to reconcile with collision-less particle, taking into account that most of the DM resides inside the galactic halo. These issues imply that in order to tackle the DM phenomenon, we need to refine the adopted DM scenario.

At this point, it is important to point-out that besides the structure of present-day DM halos, also the DM halo evolution through cosmic-time may reflect the nature of DM particle: in $\mathrm{\Lambda CDM}$ scenario, the collision-less nature of dark particle is itself imprinted in the dark halo evolution with time <cit.>, which is well followed by hydrodynamical simulations <cit.>. Also in the other scenarios of DM such as warm dark matter, Boson condensates, interacting DM, ultra light axions etc., and alternative models of DM such as Modified Newtonian Dynamics (MOND), one can find correspondence between properties of DM halos around local objects and those at high-redshift (hereafter high-$z$).

Along these lines of thought  we must remind that, until recently, the  determinations of dark matter halo density profiles was impossible: the kinematical measurements lacked of sufficient  sensitivity and spatial resolution. But, in the last decade, advanced use of integral field units (IFUs) in galaxy surveys has opened the several possibilities of studying the spatially resolved kinematics of high-redshift galaxies. For example, surveys with the Multi-Unit Spectroscopic Explorer (MUSE: ), K-band Multi-Object Spectrograph (KMOS: ), the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI: ), and Atacama Large Millimeter/submillimeter Array (ALMA: ). In particular, various resolved studies related to galaxy morphology, rotation curves, and state of inter-stellar medium has remarkably extend our view towards exciting models of galaxy evolution, for details see detailed discussion by <cit.> and <cit.>. Moreover, one recent study by <cit.> expands our view towards dark matter cusp to core transformation using nine star-forming galaxies at $z\approx 1$. In our (this) work, we focus on detailed mass modelling of rotation curves of $z\sim 1$ star-forming galaxies using a large dataset, in the context of dark matter halo and their redshift evolution.

Recalling our previous works, especially in <cit.>, we have used KMOS Redshift One Spectroscopic Survey (KROSS) data to derive the kinematics of large sample of star-forming disc-like galaxies at $z\sim 1$. We have corrected such kinematics to proper rotation curves by correcting the issues related to measurements, such as, beam-smearing[Without Adaptive Optics (AO), an IFU achieves only $0.5\arcsec -1.0 \arcsec$ spatial resolution, whereas, a galaxy from $z\sim 1$ has a typical angular size of $2\arcsec - 3\arcsec$. The finite beam size causes the line emission to smear on the adjacent pixels. This effect is referred to as `Beam Smearing', which under estimates the rotation velocity and overestimates the velocity dispersion.] and pressure-support[the kinematics of the galaxies are derived using the emission lines like H$_\alpha$ or [OIII]. These emission lines arise from the gaseous disc around the stars or the ISM. If the ISM is highly turbulent, then the emission also experiences a turbulence, i.e., radial force against gravity. This turbulence/force scales with the gas density and velocity dispersion. Since the density and velocity dispersion both decrease with increasing radius, this creates a pressure gradient, i.e., $F_{P} \propto -dP/dr$ (where $P \propto \rho \sigma^2$). The resulting radial force supports the disc and makes it rotate slower than the actual circular velocity, which might result in declining the RCs and potentially underestimate the dynamical masses] conditions. In particular, we have used 3D-forward modelling of datacubes using $^{3D}$Barolo <cit.>, which provides somewhat better constraints on kinematics than 2D-kinematic modelling approaches. For pressure support, we employed the Pressure Gradient Correction (PGC) method, which make use of full information available in datacube and therefore provides better results than pressure support correction assuming isotropic and constant velocity dispersion, for details we refer reader to <cit.>. This allows us to have a large sample of more than 200 accurate Rotation Curves (RCs) of disc-like galaxies from 6.5 Gyrs (lookback time). This large amount of data can help us open a new  portal to the Nature of Dark Matter, as well as, allows us a powerful line of investigation to test the $\Lambda$CDM collision-less scenario.

In the follow-up study, <cit.>, we employed a model independent approach of dark matter determination using the RCs established in <cit.>. Where, the  stellar mass ($M_*$) of the objects were estimated by using mass-to-light ratios derived from fitting the spectral energy distribution of the galaxies, the total gas (molecular and atomic) mass by means of  the scaling relations <cit.>, and the dynamical mass $G^{-1} V^2 \ R$  from the rotation curves, evaluated at different scale radius (disk radius: $0.59R_e$, optical radius: $1.89R_e$, and outer radius: $2.95R_e$ ). These measurements provide us with the  DM fractions ($f_{DM} = 1 - M_{bar}/M_{dyn}$) at the above radii.

In the present work, we investigate the kinematics of our same star-forming disc-like galaxies from $z\sim 1$, and dynamically mass-model their RCs. The latter is model dependent approach of deriving the physical quantities such as dark matter structural properties & its fraction, stellar masses, effective radii, and bulge masses. Thus, comparing the two studies, <cit.> and current work, has allowed us to check the robustness of our dynamical models. Through the mass-modelling techniques, we have also determined the structural properties of DM halos and checked the scaling relations between dark and baryonic matter those emerged in local star-forming disc galaxies. In the end, we compare our work with the recent study of <cit.>, in which the author dynamically mass model the galaxy kinematics, which is derived from the one-dimensional software Slit Extraction Method. The latter method is known to bias the kinematic extraction, as shown in <cit.> and <cit.>. Moreover, <cit.> corrects the pressure support with <cit.>, which underestimates the rotational velocity, as shown in <cit.>. Despite the significant discrepancies between <cit.> and our work, we compared the two studies in terms of the dark matter fraction and its structural properties.

This article is organized as follows: Section <ref> contains a brief discussion on the data employed, in Section <ref> we provide the mass models to disentangle the baryons and DM component from RCs and describe their fitting techniques. In Section <ref> we show the results and critically analyse them. In Section <ref> we compare our results with current state-of-art simulations. Finally in Section <ref> and Section <ref>, we discuss the results and provide a summary with future prospective. 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}}$.

§ DATA

We studied KROSS parent sample previously studied by <cit.> and others <cit.>, this sample contains 586 star-forming galaxies (SFGs) from  $z\sim 1$.  We have selected 344 objects out of 586, on the basis of integrated $H\alpha$ flux cut ($F_{H\alpha}>2\times 10^ {-17} \  \mathrm{[erg \ s^ {-1} \ cm^ {-2}]} $) and inclination angle ($25^{\circ} \leq \theta_i \leq 75^{\circ} $). The chosen flux and inclination cuts ensure the sufficient signal-to-noise data (S/N) and reduces the impact of extinction[In the highly inclined system ($\theta_i>75^{\circ}$) observed flux extinct due to extinction, which suppress the rotation velocity upto a few times $R_e$ <cit.>. On the other hand, in face-on galaxies ($\theta_i<25^{\circ}$) the rotation signal drops below the observational uncertainties. Therefore, to be conservative we down-select the sample for $25^{\circ} \leq \theta_i \leq 75^{\circ} $.] during the kinematic modelling.  The intrinsic characteristic of the selected sample is the following (given with respect to TableA1 of H17:): 1)AGN-flag is zero i.e., no evidence for an AGN contribution to the $H_\alpha$ emission-line profile; 2)H17 Quality-flag 1, 2, and 3,  i.e., only $H_\alpha$-emission line detected objects ($S/N>3$). We adopted the values of requisite quantities, such as effective radii ($R_e$). stellar mass ($M_*$), and redshift ($z$), from <cit.> and <cit.>. Moreover, the reddening corrected $H_\alpha$ star-formation rate ($SFR_{H\alpha}$) of the sample are adopted from <cit.>. We remark that the parent sample is selected in such a way that it does not preferentially contain galaxies in a merging or interacting state, and include good representatives of main sequence SFGs at $z\sim 1$ <cit.>.

§.§ Individual RCs

For extracting kinematics of selected sample (344 SFGs), we exploit publicly available datacubes (data availability[KROSS: http://astro.dur.ac.uk/KROSS/data.html]), and kinematics of these objects were re-derived from KROSS datacubes employing 3D-BAROLO code <cit.>. The latter takes into account the beam smearing correction in 3D-space and provides the moment maps, the stellar surface brightness profile, the rotation curve (RC), the dispersion curve (DC) along with the kinematic models. After running the 3D-Barolo, we assess the quality of the objects, which is defined by taking into account the signal-to-noise (S/N), and limitation of 3D-Barolo (e.g. inclination). The best quality objects, i.e., Quality-1 & 2,  are the one that have $\mathrm{S/N}\geq3$ and inclination $30^{\circ} \leq \theta_i \leq 75^{\circ}$.  Apart from this, our RCs were corrected for pressure gradients, which likely affect the kinematics of high-$z$ galaxies. Out of 344 KROSS objects, we choose to analyse only 256 rotation dominated Quality-1 & 2 objects, referred to as Q12 sample, for details we refer the reader to <cit.>. In the following work, we have used this Q12 sample and their RCs to analyse the dark matter distribution.

§.§ Co-added RCs

In this work we will investigate the individual RCs of Q12 sample, as well as their co-added RCs (or stacked RCs). The process of co-addition, that we briefly recall here, is accurately described in <cit.> and <cit.>. That is, the individual RCs are co-added according to their circular velocities $V_c$ (computed at $R_{out}=2.95 \ R_e$), keeping 15 RCs per bin, except for the highest velocity bin. Then, we bin the galaxies radially per $2.5 \ kpc$ corresponding to the binning scale of 3D-Barolo, see <cit.>. For the binning, we have used standard weighted mean statistic given by: $\bar{X} = \frac{\Sigma_{i=1}^{n}x_i \times w_i}{\Sigma_{i=1}^{n}w_i}$, where, $ x_i=data; \ w_i=1/error^2; \ \bar{X}=binned \ data$. The errors on binned data are root mean square error, computed as $\delta_{v}^{r_{i}} = \sqrt{\sum([(\delta_{v_{i}}^r)^2 + (\sigma_{v}^{r_{i}})^2])/N}$, where $\delta_{v_{i}}^r$ is the individual error on the velocities per radial bin and $\sigma_v^{r_{i}}$ scatter per radial bin. This process yields 16 coadded RCs, each one have minimum 5 independent data-points, see Figure <ref>, and for more details on the binning process we refer the reader to <cit.>. While binning the RCs, we also co-add & bin the physical quantities (e.g., $M_{\mathrm{D}}, R_{\mathrm{D}}$) associated with individual objects per bin using RMS statistics and the errors are standard errors. These binned parameters are used in RC modelling either as a initial guess (when the parameters have a flat prior) or as a central value in the case of a Gaussian prior. In the following text, the binned parameters are represented by a tilde character (e.g. $\tilde{M}_{\mathrm{D}}$) and co-added plus binned RCs are referred to as co-added RCs.

We remark that the main advantage of the co-added RC study is that it allows the statistical investigation of a particular class of objects, especially when the data quality is low and the number of observations are limited. The latter is currently the case for  galaxies. Moreover, $H_\alpha$ emission is the only way to trace the kinematics of  galaxies in the optical (and near-infrared) wavelength range. At the same time we know, $H_\alpha$ originates from star-forming regions, i.e., the interstellar medium, which is turbulent at . Therefore, the RCs are generally perturbed (or wiggled ). In these mixed situations, having poor data quality, and limited possibilities, the coaddition of RCs acts as a rescue tool, providing a smooth circular velocity of the galaxies, and ignoring the random fluctuations in the data. However, the coaddition of RCs also has some drawbacks, for example <cit.> shows the diversity in the RCs of local dwarf disc galaxies, that can not be evident in the case of co-added RCs. Despite the latter aspect, and taking into account the situation of  observation, we propose to study co-added RCs, as well as individuals.

Furthermore, before performing the coaddition of the RCs, we visually inspect each individual RC, which shows that the individual RCs of each velocity bin have more or less the same radial profile. Quantitatively, the data points of each velocity bin have (1) a maximum variance of $\sim0.02$ dex, and (2) a maximum dispersion of $\sim 0.1$ dex. This is very similar to the study performed by <cit.> for local SFGs. Thus, our method of co-addition is applicable and certainly necessary for  data.

§ METHODOLOGY

The rotation curve (RC) of a galaxy gives us the total circular velocity of the system as a function of its radius. This $V_c(R)$ is composed of the contribution of baryons (stars and gas) and dark matter. Therefore, $V_c(R)$ can be described as a sum in quadrature of:

\begin{equation}
\label{eq:vtot}
V^2_c(R) = V^2_{\mathrm{D}}(R) + V^2_{\mathrm{bulge}}(R) + V^2_{\mathrm{H2}}(R) + V^2_{\mathrm{HI}}(R) + V^2_{\mathrm{DM}}(R) ,
\end{equation}

were $V_{\mathrm{D}}(R), V_{\mathrm{bulge}}(R) , V_{\mathrm{H2}}(R)$ and $V_{\mathrm{HI}}(R)$ are the circular velocity profiles of the stellar disc, stellar bulge, molecular (H2) and atomic (HI) gas discs, together providing the baryon's contribution to the circular velocity and $V_{\mathrm{DM}}(R)$ gives the dark matter contribution. In the following sections, we present an approach to disentangle the aforementioned contributions from observed RCs. Note that we have corrected the observed RCs from inclination and pressure support, i.e. RCs are very close to intrinsic circular velocity curves.

§.§ Mass Models

Here we explain the various components of our mass modelling.

§.§.§ Stellar and Gas Disc
The stellar and gaseous components of local star-forming discs, dwarf-discs and low surface brightness galaxies are well described by <cit.> exponential disc model <cit.>. Here, we assume that the distribution of stars and gas in $z\sim 1$ star-forming disc galaxies are also in exponential discs, so that the surface densities:

\begin{equation}
\label{eq:freeman-density}
\Sigma_{\_}(R) \propto \frac{M_{\_}}{R_{scale}} \exp(-R/R_{scale}),
\end{equation}

where, $M_{\_}$ and $R_{scale}$ are the total mass and the scale length of the different components (stars, H2, HI), respectively. Note that stellar mass in the disc (without bulge) is denoted by $M_{\mathrm{D}}$, while the contribution of the total stellar mass is represented by $M_{\mathrm{star}} = M_{\mathrm{D}} + M_{\mathrm{bulge}}$. Given the density distribution of the stars and the gas, their contribution to the circular velocity of the disc can be expressed as follows:

\begin{equation}
\label{eq:Vd}
V_{\_}^2(R)= \frac{1}{2} \Big( \frac{GM_{\_}}{R_{scale}} \Big) \ (x^2) \  [I_0K_0 - I_1K_1],
\end{equation}

where, $x = R/R_{scale}$ and $I_n$ and $K_n$ are modified Bessel functions computed at $1.6x$ for stars and $0.53x$ for gas <cit.>.

§.§.§ Bulge Component
First, we would like to emphasise that our RCs are not well resolved within the 5 kpc where we expect the bulge to contribute significantly, i.e., the bulge is unresolved in our data. Thus, we assume that the bulge is a point mass, so that its contribution to the circular velocity becomes:

\begin{equation}
\label{eq:Vb}
V_{\mathrm{bulge}}^2(R)= \frac{G\ M_{\mathrm{bulge}}}{R},
\end{equation}

where, $M_{\mathrm{bulge}}$ is total bulge mass, and can be estimated using bulge-to-total ratio $B/T = M_{\mathrm{bulge}}/M_{\mathrm{star}}$.

§.§.§ Dar Matter Halo Component

In this work, we test two widely used DM halo models, namely Burkert <cit.> and NFW <cit.> profiles.

Burkert Halo: Local spirals and low surface brightness galaxies indicate towards the existence of central DM cores, i.e. $\rho_{inner} \propto const$ <cit.>. This observed DM-profile is well fitted by the <cit.> halo, which possesses a double power law in the DM-density, i.e., at small radii $\rho \propto R^{0}$ and at larger radii $\rho \propto R^{-3}$. Such a DM-density distribution can be expressed as:

\begin{equation}
\label{eq:BurkertDensity}
\rho (R) = \frac{\rho_{_0}}{ \Big(1+ \frac{R}{r_{_0}} \Big) \Big( 1+ \frac{R^2}{r_{_0}^2}  \Big)   ,      }
\end{equation}

where $\rho_{_0}$ and $r_{_0}$ are the central DM core density and core radius, respectively. Assuming spherical symmetry, the mass profile of Burkert DM halo follows as:

\begin{equation}
\label{eq:BurkertMass}
M_{\mathrm{DM}}^{\mathrm{Burk}} (R) = 4 \pi \rho_{_0} r_{_0}^3 \Big[ \ln \Big(1+\frac{R}{r_{_0}} \Big) - \arctan \Big(\frac{R}{r_{_0}} \Big) + 0.5 \ \ln \Big(1+\frac{R^2}{r_{_0}^2} \Big)  \Big],
\end{equation}

NFW Halo: In the standard collisionless cold dark matter paradigm, the current cosmological simulations predict a cuspy DM distribution in the center, i.e., $\rho_{inner} \propto R^{-1}$. This type of DM profile is well approximated by NFW halo, which is again a double power law but it has $\rho \propto R^{-1}$ at small radii and $\rho \propto R^{-3}$ at larger radii. Such a DM-density distribution can be expressed as:

\begin{equation}
\label{eq:NFWDensity}
\rho (R)  = \frac{\rho_{_s} }{\Big(\frac{R}{r_{_s}} \Big) \ \Big( 1+\frac{R}{r_{_s}} \Big)^2}
\end{equation}

where $\rho_{_s}$ and $r_{_s}$ are, respectively, the characteristic density and scale radius of the DM distribution. Assuming spherical symmetry, the mass profile of NFW DM halo follows as:

\begin{equation}
\label{eq:NFWMass}
M_{\mathrm{DM}}^{\mathrm{NFW}} (R) = 4\pi \rho_{_s} r_{_s}^3 \Big[\ln \Big(1+\frac{R}{r_{_s}} \Big) - \frac{\frac{R}{r_{_s}}}{1+ \frac{R}{r_{_s}}}   \Big],
\end{equation}

The circular velocity of the above DM profiles can be written as:

\begin{equation}
\label{eq:Vhalo}
V_{\mathrm{DM}}^2 (R) =  G \frac {\ M_{\mathrm{DM}} (<R)}{R},
\end{equation}

where $M_{\mathrm{DM}} (<R)$ is enclosed DM mass within radius $R$ given in Equation <ref> & <ref>. Above formalism allow us to compute the various structural properties of DM halo. For the purpose of this work, we define virial radius ($R_{\mathrm{vir}}$), virial mass ($M_{\mathrm{vir}}$), and concentration ($c_{vir}$) parameter of halo. The virial radius is the distance within which the average halo density is 200 times the critical density ($\rho_{\mathrm{crit}} = 9.3 \times 10^{-30} gm \ cm^{-3}$) of the universe, it can easily be computed numerically. The virial mass is defined as:

\begin{equation}
\label{eq:Mvir}
M_{\mathrm{vir}} = \frac{4}{3} \pi \ R_{\mathrm{vir}}^3  \  200 \  \rho_{\mathrm{crit}},
\end{equation}

In case of NFW halo, the ratio of virial radius to scale radius gives the halo concentration, $c_{vir}=R_{\mathrm{vir}}/r_{_s}$. Furthermore, once the various components of RC are disentangled, then, given the information on $V_{\mathrm{DM}}(R)$ and $V_{c}(R)$, one can also estimate:

\begin{equation}
\label{eq:fdm}
\fdm (<R) = \frac{V_{\mathrm{DM}}^2(R)}{V_{c}^2(R)},
\end{equation}

which gives us the DM fraction within the radius $R$.

Coadded and binned Rotation Curves, the so-called stacked RCs.  The individual and binned RCs are shown by blue thin curves and hexagon connected with the dark blue line, respectively. The central velocity ($\tilde{V}_c$) and the number of data points per bin are printed in the legend of each panel. For example, the notation $\tilde{V}_c:82$ (15; 64) means that the central circular velocity of the bin is 82 $\mathrm{km \ s^{-1}}$, contains 15 RCs and a total of 64 data-points. Note that here $\tilde{V}_c$ is computed at $R_{out}$ ($=2.95 \ R_e$). Note that here $V_c$ is the circular velocity computed at $R_{out}$.

Individual RC mass modelling. Row-1: Posterior distributions (MCMC output) of the estimated parameters: $r_{_{0/s}}$ and $\rho_{_{0/s}}$, for Burkert (red) and NFW (blue) halo, left and right, respectively. The vertical and horizontal gray dashed lines show the initial guess of parameters: $r_{_{0/s}}$ and $\rho_{_{0/s}}$. The full posterior distribution of parameter space $\vec{\theta}$ of these example cases are shown in Figure <ref>. Row-2 & 3 shows the best fit to the RCs and their velocity decomposition for Burkert and NFW halos, respectively. In RC panels, the blue square with error bars represents the observed data, the blue solid line shows the best fit to the data, and the blue shaded area represents the $3\sigma$ error in the fit. The decomposed velocity of stars, molecular gas, atomic gas, bulge, and dark matter is represented by yellow, violet, light blue, brown, and black dotted dashed lines, respectively. The vertical orange line shows the effective radius and grey shaded area represents the unresolved region of RC. The best-fit value of parameter space $\vec{\theta}$ are given in the legend of the figures. All the parameters are printed in log, with physical units as following: $M_{\mathrm{D}} \ [\mathrm{M_\odot}]$, $R_{\mathrm{D}} \ [\mathrm{kpc}]$, $M_{\mathrm{bulge}} \ [\mathrm{M_\odot}]$, $r_{_{0/s}} \ [\mathrm{kpc}]$ and $\rho_{_{0/s}} \ [\mathrm{gm \ cm^{-3}}]$.

Co-added RC mass modelling. Row-1: Posterior distributions (MCMC output) of the estimated parameters: $r_{_{0/s}}$, $\rho_{_{0/s}}$ and $M_{bulge}$, for Burkert (red) and NFW (blue) halo. The vertical and horizontal gray dashed lines show the initial guess of parameters: $r_{_{0/s}}$, $\rho_{_{0/s}}$ and $M_{bulge}$. The full posterior distribution of parameter space $\vec{\theta}$ of these example cases are shown in Figure <ref>. Row-2 & 3 shows the best fit to the RCs and their velocity decomposition for Burkert and NFW halos, respectively. In RC panels, the blue square with error bars represents the observed data, the blue solid line shows the best fit to the data, and the blue shaded area represents the $3\sigma$ error in the fit. The decomposed velocity of stars, molecular gas, atomic gas, bulge, and dark matter is represented by yellow, violet, light blue, brown, and black dotted dashed lines, respectively. The vertical orange line shows the effective radius and grey shaded area represents the unresolved region of RC. The best-fit value of parameter space $\vec{\theta}$ are given in the legend of the figures. All the parameters are printed in log, with physical units as following: $M_{\mathrm{D}} \ [\mathrm{M_\odot}]$, $R_{\mathrm{D}} \ [\mathrm{kpc}]$, $M_{\mathrm{bulge}} \ [\mathrm{M_\odot}]$, $r_{_{0/s}} \ [\mathrm{kpc}]$ and $\rho_{_{0/s}} \ [\mathrm{gm \ cm^{-3}}]$.

Parameter Range      Case-1 (Individual RCs)     Case-2 (Co-added RCs)

(log scale, Units)     Prior: central value, $\sigma$ (dex)     Prior: central value, $\sigma$ (dex)

$8.0 \leq \log( M_{\mathrm{D}} \ \ [\mathrm{M_\odot}]) \leq 12.5$     Gaussian Prior: $0.9 M_{\mathrm{star}}$, $\sigma_{M_\mathrm{star}}$      Gaussian Prior: $0.9 \tilde{M}_{\mathrm{star}}$, $\tilde{\sigma_{M_{\mathrm{star}}}}$

$-1.7 \leq \log(R_{\mathrm{D}} \ [\mathrm{kpc}]) \leq 1.7$     Gaussian Prior: $R_D$, $\sigma_{R_D}$      Gaussian Prior: $\tilde{R}_D$, $\tilde{\sigma_{R_D}}$

$5.0 \leq \log(M_{\mathrm{bulge}} \ [\mathrm{M_\odot}]) \leq 11.5$     Gaussian Prior: $0.1 M_{\mathrm{star}}$, $\sigma_{M_\mathrm{star}}$       Flat Prior

$-2.0 \leq \log(r_{_{0/s}} \ [\mathrm{kpc}]) \leq 2.0$     Flat Prior     Flat Prior

$-26 < \log(\rho_{_{0/s}} \ [\mathrm{gm \ cm^{-3}}]) < -18$     Flat Prior     Flat Prior

Mass modelling parameter range and their prior details. Note that Gaussian prior keep true photometric uncertainty (or dispersion $\sigma$) around the central value, where central values are computed from the high-resolution photometric data <cit.>. In case of co-added RCs, quantities are represented by tilde over-head (e.g. $\tilde{R}_D$).

§.§ Rotation Curve Modelling

To perform the mass modelling of RCs, we define a $\chi^2$ test-statistic on our observed kinematics $V_o(R)$ and modelled kinematics $V_m(R)$ as:

\begin{equation}
\label{eq:loglikehood}
\mathcal{L}_{\textrm{kin}}  \equiv  \prod_{k = 1}^{\textrm{N}} \frac{1}{\sqrt{2\pi} \ \Delta V_o(R_{_k})} \exp \Big [ -\frac{1}{2} \Big( \frac{V_o(R_{_k}) - V_m (R_{_k})}{\Delta V_o(R_{_k})} \Big)^2  \Big],
\end{equation}

For each bin $k\leq \mathrm {N}$, we compare the theoretical prediction and observed dataset by including the observational uncertainty $\Delta V_o(R_{_k})$ at each step. Our fitting procedure is based on Bayes' theorem:

\begin{equation}
\label{eq:bayes}
\mathcal{P}\left(\textrm{Model}(\vec{\theta}) \ | \ \textrm{data} \right) \propto \mathcal{L}_{\textrm{kin}}\left( \textrm{data} \ | \ \textrm{Model}(\vec{\theta}) \, \right) \ \mathcal{P}_{0}\left(\textrm{Model}(\vec{\theta}) \, \right),
\end{equation}

where the posterior probability density function (p.d.f.) is sampled from the product of the prior probability distribution assigned to the set of model parameters $ \vec{\theta}$, with the  likelihood function reported in Equation <ref>, up to the overall normalization defining the so-called  evidence of the model. The general model of RC under the Bayesian inference is defined by five parameters $\vec{\theta} = \{M_{\mathrm{D}}, R_{\mathrm{D}}, M_{bulge}, r_{_{0/s}}, \rho_{_{0/s}} \}$. Here, three parameters namely $ M_{\mathrm{D}}, R_{\mathrm{D}}, M_{bulge},$ constrain the baryons, and two parameters namely $r_{_{0/s}}, \rho_{_{0/s}}$ describe the DM halo[In $\rho_{_{0/s}}$ and $r_{_{0/s}}$, subscript $0$ gives the core radius and core density in case of Burkert halo, while subscript $s$ defines the scale radius and characteristic density of NFW halo. ]. The details of these parameters in various cases (studied in this work) of mass-modelling is given in Table <ref>. By employing Equation <ref> & <ref>, we perform the Markov Chain Monte Carlo (MCMC) analysis exploiting RC dataset discussed in Section <ref>. To perform the MCMC we make use of emcee package <cit.>. For each object, we allow 100 walkers to evolve for 10,000 steps, which gives us $10^6$ samples for analysis. We remove the first half of the sample to account for the burn-in period and the remaining are used in drawing the posteriors. Let us remark, to satisfy the Equation <ref>, in addition to the five fitting parameters given by $\vec{\theta}$, we also need the gas mass ($M_{ HI }$ and $M_{H2}$) and its scale length ($R_{ HI }$ and $R_{H2}$). Here, $M_{ HI }$ and $M_{H2}$ are derived from scaling relation <cit.> using the stellar mass computed in each MCMC run (for fixed SFR and $z$), i.e. gas mass depends on best-fit stellar mass. The gaseous disk lengths, $R_{ HI }$ and $R_{H2}$, are fixed in the fitting procedure and determined beforehand in <cit.>.

§.§ Fitting Details

As we briefly mentioned above, in fitting procedure we have five free parameters, given by $\vec{\theta} = \{M_{\mathrm{D}}, R_{\mathrm{D}}, M_{bulge}, r_{_{0/s}}, \rho_{_{0/s}} \}$, where three parameters, namely $ M_{\mathrm{D}}, R_{\mathrm{D}}, M_{\mathrm{bulge}},$ constrain the baryons and two parameters, namely $r_{_{0/s}}, \rho_{_{0/s}}$ describe the DM in case of Burkert and NFW halo. The exploration range of each parameter is motivated by mass modelling of local galaxies <cit.>, given in Table <ref>. As mentioned above, we perform RC modelling in two cases: Case-1) individual RCs, Case-2) co-added RCs. In both cases, initial guess on parameters is same and as following $\vec{\theta} = \{10, 0.3, 8.5, 1, -23.5\}$ (in dex). The prior on the dark matter parameters are kept flat in both cases. On the other hand, the baryonic parameters have different conditions on prior in different cases, each motivated by the need to treat the data robustly and without biases. Particularly, in Case-1 we allow the Gaussian prior for $M_{\mathrm{D}}, R_{\mathrm{D}}$, and $M_{bulge}$ that keep Gaussian dispersion ($\sigma$) of $\sim 65\%$ around the central value. Where, the central values for $M_{\mathrm{D}}=0.9M_{\mathrm{star}}$, which at first approximation excludes the bulge. The value of $M_{\mathrm{star}}$ and $R_{\mathrm{D}}$ are taken from <cit.>, and the bulge mass $M_{bulge}$ is $1/10^{th}$ of the stellar mass, motivated by local SFGs <cit.>. In Case-2, we allow the Gaussian prior only on $\tilde{R}_{\mathrm{D}}$, $\tilde{M}_{\mathrm{D}}$ having dispersion of $\sim 65\%$ around the central value, and all other parameters keep a flat prior. In this case central values are the quantities binned in the co-addition process of RCs, as discussed in Section <ref>. Similar to the Case-1, central value of $\tilde{M}_{\mathrm{D}}=0.9 \tilde{M}_{\mathrm{star}}$.

Note that our modelling approach also requires other parameters in the fitting process, namely $M_{\mathrm{H2}}$, $M_{\mathrm{HI}}$, $R_{\mathrm{H2}}$, and $R_{\mathrm{HI}}$ (see Equation <ref>). In particular, $M_{\mathrm{H2}}$ and $M_{\mathrm{HI}}$ are calculated in each MCMC run for a given star formation rate ($M_{H2, HI }(SFR, M_{\mathrm{D}}, z)$) using the scaling relations given by <cit.>, respectively. The molecular gas scale length $R_{\mathrm{H2}}$ is determined from the $H\alpha$ surface brightness distribution, since ionized gas around a star is surrounded by molecular gas clouds. For a detailed calculation of $M_{\mathrm{H2}}$, $M_{\mathrm{HI}}$ and $R_{\mathrm{H2}}$, we refer the reader to <cit.>. Moreover, studies of local disc galaxies have shown that the surface brightness of the HI disc is much more extended than that of the H2 disc <cit.> see also <cit.>. Therefore, we assume $R_{\mathrm{HI}}=2 \times R_{\mathrm{H2}}$, which is a rough estimate and not effecting DM fraction. Note that these additional parameters are fixed in the modelling procedure to the best of our knowledge.

Furthermore, we perform the mass modelling on individual RCs, that we have referred to as Case-1. We find that there are 28 perturbed RCs that we could not model, i.e., we mass-model 228 individual RCs. After the mass modelling of these RCs, we visually inspected each MCMC output (posteriors) and found that the 92 objects have hit the prior boundaries, either for DM parameters or for baryons (in particular $M_{\mathrm{D}}$). This tells us that the fitting of these 92 objects is prior-driven, so we have discarded them from the analysis. Thus, we end up analysing only 136 individual RCs. To verify Case-1, in which we discard substantial amount of RC-fits, we also mass-modelled the co-added RCs, which we have referred to as Case-2. In Case-2, we use up all the 228 RCs, making 16 co-added RCs (see Section <ref>), which are mass modelled using the aforementioned prescription.

§ RESULTS & ANALYSIS

As described in Section <ref> & <ref>, by inputting the individual RCs (Case-1) and the co-added RCs (Case-2), we constrain the parameter space $\vec{\theta}$ and disentangle the various  components of each RC. With this information, we have estimated the quantities $R_{\mathrm{vir}}$, $M_{\mathrm{vir}}$, $c_{vir}$, and $\fdm$. In the following section, we discuss these results.

§.§ Mass Modelled RCs

Here, we present the results of mass modelling in Case-1 (individual RCs) and Case-2 (co-added RCs). Some examples of their MCMC outputs (posteriors) and best-fit to the RCs are shown in Figure <ref> & <ref>, respectively. We note that in both cases the Burkert and NFW halo models provide a good and similar fit to the data. The rationale for this is that both halo models predict a similar DM density profile in the outskirts of the galaxies and differ only in the very inner region ($<2$ kpc), where our data are unresolved. Therefore, we cannot exclude one model over the other by just looking at the rotation curve fitting.

Next, we focused on comparing the physical properties of galaxies, such as total stellar and bulge masses. In particular, we compare the physical properties derived from Case-2 because this case gives us a smooth distribution of RCs that ignores the random fluctuations (perturbations) caused by bad data points, which makes the mass modelling more impeccable than Case-1. Therefore, in the left panel Figure <ref>, we compare the dynamically derived stellar masses from Burkert and NFW halo of Case-2 with the photometric stellar masses. As we can see, the dynamical stellar masses from the Burkert halo agree with the photometric stellar masses with an intrinsic scatter of 0.3 dex. In contrast, the dynamical stellar masses derived from the NFW halo are a factor of 0.5 lower than the photometric stellar masses, with an intrinsic scatter of 0.43 dex. This suggests that the NFW halo most likely suppresses the baryon content, which could be problematic while deriving the DM fraction.

Moreover, in both cases (with both halo models) we have difficulty in constraining the bulge mass, and often end up with a very extended bulge mass distribution. However, the bulge is never heavier than $0.3 \times M_{\mathrm{star}}$ (see right panel Figure <ref>), nor does it affect the DM fraction (see Figure <ref>). In fact, DM has emerged as a dominant quantity, especially in the outskirts[Here outskirts means the region of galaxies beyond it's effective radius.], which is very similar to the local SFGs. We also find that for low-mass objects ($M_{\mathrm{D}}< 10^{9.5} \ \Ms$), the atomic gas kinematics (HI) dominate the baryonic component, while for massive objects the stellar kinematics prevail.

Left panel: Comparison of the photometric and dynamical stellar masses in Burkert (blue) and NFW (pink) halo profiles for co-added RCs. The black solid line shows the one-to-one relations. We note that the dynamical stellar masses obtained in the possession of the Burkert halo agree very well with the photometric stellar masses. In contrast, the dynamical stellar masses obtained with the NFW halo are lower than the photometric stellar masses by a factor of 0.5. Right panel: Dynamically derived bulge masses versus stellar masses, in the case of Burkert (blue) and NFW (pink) halo profiles for co-added RCs. We compare our results with a pioneering study in this field: <cit.>, and we find that in our sample $B/T$ is never higher than $\sim 0.3$, which is different from what is reported by <cit.> for the same redshift range ($z=0.6-1.2$) and similar stellar mass sample.

Dark matter fraction within $R_{\mathrm{out}}$ for the Burkert and NFW halos, left and right, respectively. The color coding in both panels are the same, individual RCs: blue filled circles, co-added RCs: blue hexagons, local SFGs: coral red circles, and high-$z$ SFGs: coral red hexagons. The size of the makers in locals SFGs includes their errors. The blue shaded area represents the model-independent DM fraction study that we previously conducted in <cit.>.

§.§ DM Fraction

In Figure <ref> we show the results of DM fraction of mass modeled RCs, in Case-1 and Case-2, for Burkert and NFW halo. Firstly, we notice that both halo models in both cases have shown that the majority ($\sim 70\%$) of galaxies are DM dominated in the outskirts. Secondly, our current results are in very good agreement with our previous study <cit.>, where we determine the DM fraction using a model-independent approach. This confirms that the majority of SFGs at $z\sim 1$ have dark matter dominated discs.

In case of Individual RCs, there are nearly 20% objects those appear to be deficient in dark matter (in both halo models). We also note that the DM fraction derived from the individual RC modelling retains a higher uncertainty in both halo models, while it is significantly reduced in the case of the co-added RCs. This could be due to the low resolution in the data (non-smooth $V(R)$), which is penalized by the large uncertainty. On the other hand, co-adding RCs statistically ignores the random fluctuations arising from the bad data points and provides us with a smooth RC, which can then be easily modelled dynamically, and hence keeps the relatively small uncertainty in the measurements. Besides the similarity in the $f_{\mathrm{DM}}$ measurements derived from co-added RCs and <cit.>, we also found that binning of individual $f_{\mathrm{DM}}$ in stellar mass range (keeping 10 points per bin) yields the DM fraction which overlaps the <cit.> results.

As we have previously pointed-out in Section <ref> that dynamical modelling in possession with NFW halo attempts to suppress stellar mass content, as shown in the left panel of Figure <ref>. Therefore the DM fraction obtained in this case is higher in comparison to Burkert halo case, see right panel of Figure <ref>. The latter argument is clearly visible from co-added RC (blue hexagons) analysis. This suggests that dynamical modelling with the NFW halo requires a precise baryon information, while the Burkert halo does not and is still able to constrain the physical properties of galaxies, for example stellar masses. This is the reason $\fdm$-$M_*$ relation obtained with Burkert halo matches very well with the $\fdm$-$M_*$ relation of the model-independent approach (blue shaded region). However, this can not be a piece of strong evidence to discard one model over the other because measurements from NFW models are still in good agreement within the 1$\sigma$ uncertainty.

Concentration vs. virial mass relation for individual and co-added RCs (or galaxies), blue circles and hexagons, respectively. The blue and red, solid and dashed lines, represent the $c-M_{\mathrm{vir}}$ relation of DM only simulations <cit.> and NIHAO suite of hydrodynamical cosmological simulations <cit.>, for $z= 0$ and $z\sim 1$, respectively.

DM core density vs. core radius relation of local and high-$z$ galaxies. The blue circles and hexagons represent the $z\sim 1$ individual and co-added RCs, respectively. The coral circles and stars represent the local star-forming and low surface brightness galaxies <cit.>, respectively. The brown dashed line is the best fit $\rho_{_0} \propto r_{_0}^{\alpha}$ to $z\sim 1$ data with slope $\alpha=-1.54$, and the orange shaded area represents the $1-3\sigma$ intrinsic scatter in the relation.

Dark matter fraction as a function of radius in various velocity bins of observed and simulated galaxies. The central value of $f_{ DM }(R)$ is shown by solid lines in each case. The color code is the same in all panels and is as follows: EAGLE : blue, TNG100: orange, TNG50: green, observations in the case of the Burkert halo: black, and NFW halo: brown. In the simulations, the $1\sigma$ error in each measurement is represented by a shaded area, while in the observations it is a dashed line following the same color code. The grey shaded area in each panel represents the unresolved region in the observations. The central value of the circular velocity of each bin is printed at the top of each panel, where $V_c^{ SIM }$ and $V_c^{ OBS }$ specifies the simulations and observations, respectively.

Structural properties of the DM halo at $z\sim 1$. Upper panel: DM core radius and core density as a function of circular velocity, upper and lower panels, respectively. The individual and co-added RCs measurements are represented by small blue circles and the hexagons, respectively. The comparison sample of local LSBs <cit.> and SFGs <cit.> are represented by the coral red stars and circles, respectively. The best-fit line (second-order polynomial) of the local SFGs and LSBs is represented by a red line, and the equation is printed in the upper left corner of each panel. As we can notice the local SFGs lie above in the $r_{_0} - V_c$ plane and below in the $\rho_{_0}-V_c$ plane with respect to the $z\sim 1$ SFGs. This is due to the fact that DM cores of $z\sim 1$ SFGs are smaller and denser than their local counterparts, indicating an evolution in the structure of DM halos over cosmic time, that is fairly captured in the observations. Lower panel: The ratio of DM core radius and density at $z\approx 0$ to $z\sim 1$ as a function of circular velocity. This plot shows the expansion factor of DM cores at $z\sim 1$ with respect to local SFGs. The color codes for data points are similar to those in the left panel. The blue solid line shows the best-fit to the data, obtained by fitting the second order polynomial function printed in both panels. The dashed blue lines followed by solid line shows the 1$\sigma$ intrinsic scatter in the relations.

§.§ Structural Properties of DM

With the information on the DM halo parameters: $r_{_{0/s}}, \ \rho_{_{0/s}}$, we derived the physical properties of the galaxies, such as $R_{\mathrm{vir}}, M_{\mathrm{vir}}$ and the concentration in the case of the NFW halo. Using the latter quantities, we compared the $c_{vir}-M_{\mathrm{vir}}$ relation of our sample with DM only simulation <cit.> and NIHAO hydrodynamical simulation <cit.>, as shown in Figure <ref>. We find that the concentration at $z\sim 1$ SFGs goes in the opposite direction, which challenges the simulations. In other words, both DM only and the hydro simulations predict that the halo concentration decreases with increasing redshift, (see red and blue dashed lines in Figure <ref>, while the observations show a higher concentration at $z\sim 1$ compared to $z\approx 0$, but still consistent within the $1\sigma$ uncertainty. We would like to notice the reader that the concentration is derived as $c=R_{\mathrm{vir}}/r_{s}$, where the determination of $R_{\mathrm{vir}}$ requires a huge extrapolation. Therefore, the concentration in this work is not as accurate as can be derived from HI RCs, where data are available over more than 50 kpc radius. Thus, although the NFW halo gives conflicting results with respect to the simulations, we cannot falsify its presence in nature.

On the other hand, the relation between core density versus core radius, associate with the Burkert profile, emerge very similar to the local SFGs, with a similar slope, as shown in Figure <ref>. Moreover, we find that high-$z$ galaxies have denser DM cores, as they lie above the local SFGs and low surface brightness galaxies (LSBs). At the same time, the size of their DM core is smaller with respect to local SFGs, but equivalent to the core of LSBs. This suggests that the structural properties (i.e., structure) of DM have evolved over cosmic time. In particular, DM cores at $z\sim 1$ are smaller by a factor of 0.3, which makes them denser by $\sim 1.5$ order of magnitude.

§ COMPARISON WITH SIMULATIONS

In this section, we compare the mass modelling results with the current state-of-the-art galaxy simulations: EAGLE, TNG100 and TNG50  . The details of their sample selection and RC data can be found in the Appendix <ref>. In particular, we compare our observed co-added RCs with the simulated co-added circular velocity curves. Briefly, in simulations we know exactly the total circular velocity of the galaxy as a function of radius ($V_{c}(R)$), and its components such as $V_{\mathrm{star}}(R)$,  $V_{\mathrm{\tiny {DM}}}(R)$, $V_{\mathrm{gas}}(R)$ etcetera. From observations we only know the $V_{c}(R)$, which we have mass modelled to disentangle the various components, see Section <ref>. That is, in both cases, observations and simulations, we know the $V_{c}(R)$ and $V_{\mathrm{DM}}(R)$. This allows us to compute DM fraction of galaxies as a function of radius using Equation <ref>. In Figure <ref> we show the results of $\fdm(R)$ in the aforementioned simulations. In particular, we show 12 velocity bins for which we have enough simulated data to bin and obtain the co-added RC, for details see Table <ref>. It is noticeable that all simulations show DM dominated galaxies from 5 kpc outwards. In other words, the simulated galaxies have 20% (minimum) to 50% DM within the effective radii, and they are DM dominated in the outskirts. These results are very similar to the observed SFGs at $z\sim 1$, reported previously in <cit.>. We remind the reader that galaxies within effective radii are known to have less DM fraction because baryon dominates the inner region of galaxies. If one has to look for DM, it should be found in the outskirts ($> 2 \times R_e$), as clearly shown in Figure <ref>.

To compare the simulations with the observations, we superimposed the $\fdm (R)$ of the observed co-added RCs in Figure <ref>, in the case of the Burkert (black) and NFW (brown) halos. Since the observations in the inner region are not resolved, we mark the inner region (0-4 kpc) as a gray shaded area and avoid drawing any conclusions there. Notice that we compare observed and simulated co-added RCs only when $V_c^{ SIM } \approx V_c^{ OBS }$, shown at the top of each panel. First, we note that $\fdm(R)$ in the Burkert halo case represents the simulations well within the $1\sigma$ uncertainty. On the other hand, the NFW halo yields a higher DM fraction up to the outer edges in all velocity bins. Second, NFW and Burkert halo results are very similar in low velocity bins ($V_c^{ SIM }< 150 \ \mathrm{km \ s^{-1}}$) and agree with the simulations. However, when we go towards the higher velocity range, the NFW halo starts to show a deviation from the simulations, while the Burkert halo stays very close to the simulations. In particular, the Burkert halo and the simulations have almost the same DM fraction between 10-15 kpc radius, which is different from the NFW halo. In short, the Burkert profile seems to represent the simulations better and more accurately.

To obtain more quantitative difference, we emphasize to notice the DM fraction between 5-10 kpc radius in Figure <ref>. We can clearly notice that NFW halo allows more DM content than Burkert halo, which we do not even see in simulations. In other words, DM fraction derived from the Burkert halo model is close to the one seen in hydrodynamical simulations for over all stellar mass (and circular velocity) range, while NFW yields a relatively high DM fraction. The latter could be a consequence of suppressing the stellar mass which we have already anticipated from the observations in Section <ref> & <ref>,  and Figure <ref> &  <ref>.

Besides this, in observations we notice to have three co-added RCs, corresponds to $M_{\mathrm{star}} \approx 9.2, 9.4, 10.0$, with low DM fraction. We notice that these RCs have a gas-dominated outer disc, which we could not find in the simulations. Moreover, we recall that the DM fraction derived from the model-independent approach in <cit.>, overlaps with the DM fraction obtained with Burkert halo model, and so the simulations. That is, these results are strengthening our previous findings, as well as confirming the robustness of the current measurements and the quality of the KROSS data.

§ DISCUSSION

To fit the individual RCs through mass modelling approach, one needs to use the informative prior on the baryons ($M_{\mathrm{D}}, M_{bulge}, R_{\mathrm{D}}$) in advance, otherwise the parameter space $\vec{\theta}$ is highly degenerate. In the case of co-added RCs, on the other hand, we could achieve a good fit only by specifying one or two baryonic parameter, which is stellar disc scale-length and stellar mass in this work. This case also allows us to adopt the flat prior on bulge mass, which is an unknown quantity (and un-resolved), so that we can determine it dynamically, see Section <ref> and Table <ref>. Remark: when we use the NFW halo model for mass modelling of RCs, it gives a good fit to the data, but suppresses the baryon content and yields high DM fractions, which we do not see when using the Burkert halo, nor in the simulations (EAGLE, TNG100, TNG50). Based on RC fitting and structural properties of DM we can not include or exclude the NFW halo model but there exists the tension in constraining the stellar masses using NFW halo, see Figure <ref>. In principle, a galaxy can not have stellar mass more than the allowed by dynamics <cit.>. Providing the various results/comparisons of NFW and Burkert halo in the case of individual and co-added RCs, we found that our current sample favours the Burkert halo model and also agrees quite well with the simulations. Therefore, we focus on discussing the Burkert halo results in what follows.

We begin by referring to our earlier study of the DM fraction <cit.>, in which we used a model-independent approach to determine the $\fdm$ from RCs. We found that our earlier finding agrees very well with our current results of the DM fraction obtained via mass modelling (see left panel Figure <ref>). Both the current and previous $\fdm$ results also matches very well with the simulations (see Figure <ref>). Our current $\fdm$ results are also consistent with the previously studied local and high-$z$ SFGs, indicated by coral red circles and hexagons, respectively, in the Figure <ref>. Given the similarities between the observed data and the simulations mentioned above, we conclude that the current state-of-the-art simulations reproduce the observed universe well. However, for small objects with low velocities and gas dominated system, there is still some work needed to be done.

In one of our previous works <cit.>, we have shown that the progenitors of present-day SFGs at $z\sim 1$ have similar rotation curves. Given the similarity in the RCs of $z\sim 1$ and $z= 0$, and acknowledging the fact that the circular velocity is an indicator of the total mass of the system, we have concluded that the total mass within the optical extent of galaxies does not evolve from $z\sim 1$ to $z= 0$, while the baryon distribution evolves. To bring our previous confirmation and the current finding on the same level, we compared our recent results on the structural properties of DM halo in the circular velocity plane, as shown in Figure <ref>. Let us remark that $r_{_0}-V_c$ and $\rho_{_0}-V_c$ relations of local galaxies are well fitted by second order polynomials, see upper panel Figure <ref>, which are given by:

\begin{equation}
\label{eq:rc-Vc-z0}
\log r^{z\sim 0}_{_0} = 1.63 \ \log V_c^{2} - 5.24 \ \log V_c + 4.76
\end{equation}

\begin{equation}
\label{eq:rho-Vc-z0}
\log \rho^{z\sim 0}_{_0} =  -1.99 \ \log V_c^{2} +6.67 \ \log V_c -1.04
\end{equation}

Therefore, to visualize the differences in the structural properties of DM at different epochs, we used Equation <ref> and <ref> as a benchmark and have plotted the ratio: $\frac{r_{_0}^{z\sim 0}}{r_{_0}^{z\sim 1}}$ and $\frac{\rho_{_0}^{z\sim 0}}{\rho_{_0}^{z\sim 1}}$, as a function of circular velocity in the lower panel of Figure <ref>. This immediately shows us how the structural properties of DM differ at $z\sim 1$ with respect to $z\approx 0$. First, we find that the DM cores at $z\sim 1$ are smaller than their local counterparts by a median factor of 0.3. Second, the density of DM cores at $z\sim 1$ is more than an order of magnitude higher than at $z\approx 0$. It is also very clear that these two ratios correlate very well with the circular velocity of the galaxies. Thus gives us the direct relation of expansion factor for core radius and density as a function of the circular velocity of a galaxy:

\begin{equation}
\label{eq:exp-rc}
\log \Big( \frac{r_{_0}^{z\sim 0}}{r_{_0}^{z\sim 1}} \Big) = 2.1 \ \log V_c^2 - 7.3 \ \log V_c + 6.6 \ \ \ \ \ \lfloor{\pm 0.6 \ \mathrm{dex}} \rceil
\end{equation}

\begin{equation}
\label{eq:exp-rhoc}
\log \Big( \frac{\rho_{_0}^{z\sim 0}}{\rho_{_0}^{z\sim 1}} \Big) = -2.7 \ \log V_c^2 + 8.2 \ \log V_c -5.8 \ \ \ \ \ \lfloor{\pm 0.9 \ \mathrm{dex}} \rceil
\end{equation}

The Equation <ref> & <ref>, are the efficacious in relating the structural properties of DM to baryonic matter. Plotting these relations in the lower panel of Figure <ref>, clearly shows us that at $z\sim 1$ objects those possessing higher circular velocity are denser in nature, i.e. more compact in size (in terms of DM cores ), relative to $z\approx 0$ SFGs. In other words, objects with higher circular velocity expand their DM cores larger over cosmic time, relative to objects with lower circular velocity. This means that the circular velocity plays a crucial role in the galaxy evolution process. We report that the intrinsic scatter in relation <ref> & <ref>, is about 0.57 and 0.94 dex, respectively. Besides this, we also found that the galaxies in our sample have a typical DM core size of 3-10 kpc, which is a factor of 2.5 smaller than those reported in <cit.>.

Moreover, the above discussion suggests that DM is most likely responding to the baryonic processes taking place in the inner region of galaxies throughout cosmic time. This also implies that early (say, $z> 2$) DM halos are likely cuspy, which may have changed over cosmic time under the influence of galactic processes (e.g., AGN and supernova feedbacks, dynamical frictions), which has been previously suggested by many authors such as <cit.> and <cit.>. Although, the difference in the size and density of the DM-core between $z\sim 1$ and $z\approx 0$ is small, it is nevertheless remarkable and could have a very strong impact on the identification of the Nature of Dark Matter . In short, it shows a clear expansion in the structure of DM halos over cosmic timescales. Based on these results, one can easily interpret that if the DM cores are denser at $z\sim 1$, then the virialization scale (radius) of these systems should also be smaller with respect to the locals. The latter argument can be further tested using the forthcoming observations of neutral hydrogen in the outskirts of galaxies <cit.>.

Furthermore, the results obtained with the Burkert halo model are consistent with the simulations. Therefore, so far we suggest that the rotation-dominated star-forming galaxies in the simulations and observations are very likely to have a DM core at the center of the halo. However, the inference of DM cores is very indirect in this study due to limited resolution in the inner region of galaxies. Therefore, to confirm the our arguments, we still require high-resolution and good signal-to-noise data in the inner region of the galaxies. That can possibly be obtained with the upcoming James Webb Space Telescope (JWST) and in the future with the Extremely Large Telescope (ELT).

§ SUMMARY & CONCLUSIONS

In this work, we have studied the structural properties of the dark matter in $z\sim 1$ main sequence star-forming disc-like galaxies. We first mass-modelled the individual and co-added rotation curves, assuming that the baryons are distributed in the <cit.> disc and dark matter is modelled with the Burkert and NFW halo profiles. Then the mass modelling results in the Burkert and NFW cases are critically examined, interpreted, and compared with current state-of-the-art galaxy simulations. We found that the results obtained using the Burkert halo fits well with the observations and also mirror the simulations, see Figure <ref>, <ref>, <ref>, <ref>, & <ref>. Therefore, in the end, we discussed and concluded the results obtained with the Burkert halo profile, we found that:

* Our mass modelling approach has proven robust to disentangle the various components of RCs, and being able to determine physical properties of baryons such as scale length, stellar and bulge mass, and DM structural properties such as core density and radius, as well as the fraction of dark matter, for details see Appendix <ref>.

* The dark matter fraction derived from the mass modelling outcomes (using  Equation <ref>) agrees very well with the results of the model-independent dark matter fraction (obtained in ). Moreover, it is also similar to those predicted by different simulations: EAGLE, TNG100 and TNG50.

* The results of the dark matter structural properties ($r_{_0}$ and $\rho_{_0}$) show that the dark matter cores at $z\sim 1$ are on an average a factor of 0.3 smaller and more than an order of magnitude denser than those of local star-forming galaxies. This suggests an expansion in the structure of DM halos over 6.5 Gyrs.

* We also report that the expansion of DM cores is a function of circular velocity. The objects with higher circular velocity expand their core more than those with low velocities. This suggests that the circular velocity plays a crucial role in deriving the galaxy evolution.

We have shown in our previous work that the progenitors of local star-forming galaxies at $z\sim 1$ are similar in terms of total mass ($V_c^2(R) \propto M_{\mathrm{tot}}$), they are dark matter dominated objects, and their baryon distribution evolves with time as a function of their circular velocity <cit.>. Now, this pilot study of dynamical mass modelling of high-$z$ RCs is showing that the dark matter distribution also changes over cosmic timescales, and it is a function of circular velocity. Therefore, we suggest that galactic processes that cause an evolution of the baryon distribution, also have an impact on the dark matter, consequently dark matter distribution changes. In other words, dark matter responds to baryonic processes (e.g., AGN, supernova feedbacks) over the time scale of a few Gyrs ($\sim 6$). If this work is confirmed by other observations, we will have the first observational evidence of gravitational potential fluctuations in the inner region of galaxies, originally theoretically predicted by <cit.>. These evidence will have a strong implication for cross-checking simulations, as well as determining the astrophysical and particle nature of dark matter.

§.§ Acknowledgements

We thank the anonymous referee for their constructive comments and suggestions, which have significantly improved the quality of the manuscript. We thank K. Oman for providing us with the hydro-dynamical simulations and discussion/comments on the results. G.S. thanks J. Freundlich for passing the NIHAO simulations` concentration results. G.S. also thanks A. Marasco, M. Petac and S. Goswami for their fruitful discussion and various comments in the entire period of this work. 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).

§.§ Data Availability

For Physical properties of galaxies such as stellar mass, gas mass, SFR, velocity etc., we refer reader to the catalogues released with  (Flat RCs catalog ) and  (DM fraction catalog ). The measurements from co-added RC mass modelling are available in form of tables in Appendix <ref> and measurements from individual RCs are attached in the external appendix, which we will be released with the publication of this paper.

§ SIMULATIONS

Comparison of physical parameters of simulations (EAGLE: blue dots; TNG100: orange stars; TNG50: green triangles) with observations (brown filled circles). Left panel: Star formation rate as a function of stellar mass, so-called main sequence of SFGs. Middle panel: Effective radii of galaxies as a function of stellar mass, so-called mass-size relation. Right panel: Distribution of the maximum circular velocity of the simulated sample compared to $V_{out}=V_c$ of the observations.

Co-added circular velocity curves of galaxy simulations from first-to-last EAGLE, TNG100 and TNG50, respectively. The color for each panel (and plot) is given in the legend, as well as the name of the simulation and velocity bin is shown in the title of each plot.

Binned circular velocity curves of galaxy simulations from first-to-last EAGLE, TNG100 and TNG50, respectively. The color for each panel (and plot) is given in the legend, as well as the name of the simulation and velocity bin is shown in the title of each plot.

Binned RCs of all velocity bins of galaxy simulations from first-to-last EAGLE, TNG100 and TNG50, respectively. The color maps gives an idea of their increasing circular velocities.

In addition to observations, we also analyse recent state-of-the-art galaxy simulations, EAGLE <cit.> and IllustrisTNG <cit.>. For the purpose of this paper, we select central galaxies from the simulations using friends-of-friends halo finder approach <cit.> and the SUBFIND algorithm <cit.>. Gas and star particles are associated to the galaxy of their nearest neighbouring dark matter particle if belongs to one, and they are gravitationally bound to the object. The centre of a galaxy is defined as the location of its particle/cell (of any type) with the minimum gravitational potential energy evaluated over the particles belonging to the galaxy- we take this as the origin of the profiles assumed in Section <ref>.

To match the observed and simulated galaxies, we further filter the central galaxies using following selection criteria: $z= 0.85$, $M_* = 10^8-10^{11.6} \ \mathrm{M_\odot}$, SFR = $0.1-100 \ \mathrm{M_\odot yr^{-1}}$, $R_e = 0.03-30 \ \mathrm{kpc}$, $V_c = 75- 300 \ \mathrm{km \ s^{-1}} $ and  $V_{rot}/\sigma >1$. Here, $V_c$ is maximum circular velocity of the galaxies within the radius of 100 kpc radius. The $V_{rot}/\sigma$ is ratio between maximum rotation (azimuthal) velocity to the velocity dispersion. With this selection criteria, we obtained 212 EAGLE, 300 TNG100, and 551 TNG50 galaxies. We present this selected sample together with the observed data in Figure <ref>, where we show the main-sequence of SFGs, the mass-size relation, and the velocity distribution. As we can notice that our observed sample agrees very well with the simulations in terms of their effective radii, stellar masses, and circular velocities. However, the SFR of the simulated sample is lower by a factor of 0.37 dex, which is one of the well-known problems in current simulations <cit.>. In spite of this mismatch, we consider the selected simulated galaxy sample as an optimal sample for a detailed comparison with our data, as it is a representative for the population of typical star-forming, rotation-dominated galaxies at $z\sim 1$ in the $\Lambda$CDM Universe.

To compare the simulations with observations, we extract the mass profile of simulated galaxies expressed as circular velocity curve ($V_c(R) = \sqrt{GM/R}$) and its constituent components (i.e., stars, gas and dark matter). We co-add and bin the simulated circular velocity curves using the similar approach as for observations; however using the RMS statistic, which works best for simulations. Briefly, we divide our simulated galaxies into 20 velocity bins between 100-300 $\mathrm{km \ s^{-1}}$ with a step- size of 10 $\mathrm{km \ s^{-1}}$. Then each bin is co-added and binned to obtain a single circular velocity curve. We would like to note that in the simulations we know exactly the total velocity of the galaxy, as well as its different components (e.g. star: $V_*$, gas: $V_{H2}, V_{\mathrm{HI}}$, dark matter: $V_{\mathrm{DM}}$), so in the binning process we also bin these quantities. Moreover, we bin the physical quantities such as $R_{\mathrm{D}}$, $M_*$, $V_c$, etc., associated to each galaxy. In Figure <ref>, we show an example case of co-adding, binning, and the final co-added circular velocity curves from all three simulations, and in Table <ref>, we tabulate all their relevant physical quantities. These circular velocity curve data are then used in the investigation of the DM fraction, see Section <ref>. We remark, TNG100 and TNG50 simulations are same and only differs in terms of resolution. However, we still notice a little dis-similarity in the shape of their circular velocity curves. Nevertheless, in this work, we avoid making any interpretations on the circular velocity curves derived from the simulations, as this is beyond the scope of this work and will be covered in our future work, where we plan to make the synthetic datacubes of simulated galaxies, and determining the RCs using the 3D-Barolo, as suggested by <cit.> and <cit.>.

§ TESTING MASS MODELLING APPROACH

An example of synthetic co-added RCs. Top panel: an observed co-added and binned RC, which is shown as a reference. Bottom panel: co-added and binned Synthetic RCs in case of Burkert and NFW dark matter halo, left and right, respectively. The color code in all the panels are the same and given as follows: the brown square represents the individual RCs and the blue stars represents the coadded RC.

An example of MCMC posteriors of mass-modelled synthetic coadded RCs, in case of Burkert and NFW halo, left and right, respectively. Here we show the posteriors corresponds to the synthetic co-added RC shown in Figure <ref>. The best-fit parameters are printed above each posterior with $1\sigma$ uncertainty. The vertical black dashed lines show the best fit central value along with $1\sigma$ upper and lower limit, and the true values are shown by blue crosses.

To test the performance of our mass modelling approach, we verified it on synthetic RCs. In this section, we explain this procedure in detail.

§.§ Synthetic RCs

To generate the synthetic RCs, we take the reverse approach, i.e., we use the best-fitting parameters ($\vec{\theta} = [M_D, R_D, M_{\mathrm{bulge}}, r_0, \rho_0]$, gass disc parameters = [$M_{H2}, M_{ HI }, R_{H2}, R_{ HI }$], SFR, z=0.85) obtained from the mass modelling of the individual RCs. Then, the synthetic RCs are generated using the baryon and dark matter models defined in Section <ref>, and at the end we assign some noise to each RC. In this way, we have obtained 228 individual synthetic RCs, for which we know exactly the baryons and dark matter content. These synthetic RCs can thus be used to test our mass modelling approach. An example of such individual RCs is shown in the bottom panel of Figure <ref>. Note that we generate synthetic RCs for, both, the Burkert and NFW halo profiles.

§.§ Synthetic Co-added RCs

Co-addition of synthetic RCs is performed using the same strategy applied on observed RCs, defined in Section <ref>. That is, each synthetic co-added RC contain 15 synthetic individual RCs, those are binned using weighted mean statistics. During binning process, we also bin the parameters space $\vec{\theta}$ and gas disc parameters, yielding $\tilde{M}_D, \ \tilde{R}_D, \ \tilde{M}_{\mathrm{bulge}}, \ \tilde{r}_0, \ \tilde{\rho}_0, \  \tilde{M}_{H2}, \ \tilde{M}_{ HI }, \ \tilde{R}_{H2}, \ \tilde{R}_{ HI }$ for synthetic co-added RCs, which we use as a central values in mass modelling (if required). An example of synthetic co-added RC is shown in Figure <ref>.

§.§ Mass Modelling of Synthetic Co-added RCs

For the mass modelling of synthetic co-added RCs, we use the same technique as applied on observed co-added RCs, and described in Section <ref>. We examine the parameter space ($\vec{\theta}$) using the MCMC sampler. In fitting process we keep Gaussian prior for $R_D$ and $M_D$ with a 65% spread around the central values ($\tilde{M}_D, \ \tilde{R}_D,$). The rest of the $\vec{\theta}$- space keep flat prior, ranges are given in the Table <ref>. The gas mass (atomic and molecular) depends on the scaling relations <cit.> and derived in each MCMC run using the running value of the stellar mass, the fixed SFR, and $z=0.85$. That is, the best-fit gas mass depends on the best-fit stellar mass, SFR, and $z$. The lengths of the gaseous discs are fixed.

Results: In Figure <ref> lower panel, we show an example of synthetic individual and co-added RCs for bin number five (bin_5) in case of Burkert and NFW halo profiles, left and right, respectively. These synthetic RCs are shown in comparison with observed co-added RCs (upper panel) of same velocity bin. Required information related to each co-added RC is printed in the title of each RC plot. We allow the mass modelling of this (and other 15) RCs using the  aforementioned procedure. In Figure <ref>, we show the posterior of bin_5, as one can notice, our mass modelling is capable of retrieving the true $\vec{\theta}$-space within $1\sigma$ uncertainty. Therefore, we conclude that– our mass modelling approach has proven robust to disentangle the various components of RCs.

§ TABLES

In Table <ref>, we provide the physical parameters of observed galaxies (in case of co-added RCs). In particular, we provide relevant details on 1) co-added RCs (binning), 2) averaged physical quantities ($\tilde{R}_{\mathrm{D}}$, $\tilde{V}_c$, $\tilde{M}_*$, $\mathrm{\tilde{SFR}}$ and $\tilde{R}_{\mathrm{H2}}$) per bin and 3) best-fit modelled parameters space $\vec{\theta}$ and DM fraction in case of Burkert and NFW halo. In Table <ref>, we provide the relevant physical parameters ($\tilde{R}_{\mathrm{D}}$, $\tilde{V}_c$, $\tilde{M}_*$, $\tilde{R}_{gas}$, and $\tilde{f}_{\mathrm{\tiny{DM}}}$) of galaxy simulations: EAGLE, TNG100 and TNG50.

§ EXTRA FIGURES

In Figures <ref> & <ref>, we provide posterior  distribution of full parameter space  $\vec{\theta}$ for individual and co-added RCs, respectively.

col sep=comma,
string type,
every head row/.style=
before row=
19|c|Physical Parameters: OBSERVATIONS

2cBinning Details     5|c||Photometric measurements     5|c|$\vec{\theta}$: Burkert Halo     5|c|$\vec{\theta}$: NFW Halo     2|cDM Fraction

after row=,
every last row/.style=after row=,
columns/Bin/.style=column name=No., column type=l,
columns/RC_pts/.style=column name=Dpts, column type=l|,
columns/Rd/.style=column name=$\tilde{R}_{\mathrm{D}}$, column type=l,
columns/Vc/.style=column name=$\tilde{V}_c$, column type=l,
columns/Md/.style=column name=$\tilde{M}_*$, column type=l,
columns/SFR/.style=column name=$\mathrm{\tilde{SFR}}$, column type=l,
columns/RH2/.style=column name=$\tilde{R}_{\mathrm{H2}}$, column type=l|,
columns/Brd/.style=column name=$\tilde{R}_{\mathrm{D}}$, column type=l,
columns/Bmd/.style=column name=$\tilde{M}_{\mathrm{D}}$, column type=l,
columns/Bmbu/.style=column name=$\tilde{M}_{bul}$, column type=l,
columns/Brc/.style=column name=$\tilde{r}_0$, column type=l,
columns/Brho0/.style=column name=$\tilde{\rho}_0$, column type=l|,
columns/Nrd/.style=column name=$\tilde{R}_{\mathrm{D}}$, column type=l,
columns/Nmd/.style=column name=$\tilde{M}_{\mathrm{D}}$, column type=l,
columns/Nmbu/.style=column name=$\tilde{M}_{bul}$, column type=l,
columns/Nrc/.style=column name=$\tilde{r}_0$, column type=l,
columns/Nrho0/.style=column name=$\tilde{\rho}_0$, column type=l|,
columns/Bfdm/.style=column name=$f_{\mathrm{\tiny{DM}}}^{\mathrm{\tiny{BUR}}}$, column type=l,
columns/Nfdm/.style=column name=$f_{\mathrm{\tiny{DM}}}^{\mathrm{\tiny{NFW}}}$, column type=l,

Structural (and physical) parameters of observed galaxies derived from their co-added RCs. col:1-2 Bin number and total number of RCs, respectively. For example, the notation 15(67) tells that we have coadded 15 RCs those contain 67 data-points. col:3-7 Photometric measurements of physical quantities $\tilde{R}_{\mathrm{D}}$, $\tilde{V}_c$, $\tilde{M}_{\mathrm{D}}$, $\mathrm{\tilde{SFR}}$ and $\tilde{R}_{\mathrm{H2}}$, respectively. The tilde $\tilde{•}$ at the top of each quantity represents that it is a statistically measured quantity. col:8-12 Modelled physical parameters ($\tilde{R}_{\mathrm{D}}$, $\tilde{M}_{\mathrm{D}}$, $\tilde{M}_{\mathrm{bul}}$, $\tilde{r}_{_0}$ and $\tilde{\rho}_{_0}$, respectively) for Baryons and DM in case Burkert Halo. col:13-17 Modelled physical parameters ($\tilde{R}_{\mathrm{D}}$, $\tilde{M}_{\mathrm{D}}$, $\tilde{M}_{\mathrm{bul}}$, $\tilde{r}_{_s}$ and $\tilde{\rho}_{_s}$, respectively) for Baryons and DM in case NFW Halo. col:18-19 DM fraction in case of Burkert and NFW halo, respectively. Here we provide each quantity in $\log$ and their physical units are as following: $\tilde{R}_{\mathrm{D}} \ [\mathrm{kpc}]$, $\tilde{V}_c \ [\mathrm{km \ s^{-1}}]$, $\tilde{M}_{\mathrm{D}} \ [\mathrm{M_\odot}]$, $\mathrm{\tilde{SFR}} \ [\mathrm{M_\odot \ yr^{-1}}]$, $\tilde{R}_{\mathrm{H2}} \ [\mathrm{kpc}]$, $\tilde{M}_{\mathrm{bul}} \ [\mathrm{M_\odot}]$, $\tilde{r}_{_0} \ [\mathrm{kpc}]$ and $\tilde{\rho}_{_0} \ [\mathrm{M_\odot \ kpc^{-3}}]$.

col sep=comma,
string type,
every head row/.style=
before row=
19|c|Physical Parameters: SIMULATIONS

1c|Bin    3c|Number of RCs     3|c|$\tilde{R}_{\mathrm{D}}$     3|c|$\tilde{V}_c$     3|c|$\tilde{M}_*$     3|c|$\tilde{R}_{gas}$    3|c$\tilde{f}_{\mathrm{\tiny{DM}}}$

after row=,
every last row/.style=after row=,
columns/bin/.style=column name=No., column type=l|,
columns/EN/.style=column name=EAG, column type=l,
columns/T100N/.style=column name=T100, column type=l,
columns/T50N/.style=column name=T50, column type=l|,
columns/Erd/.style=column name=EAG, column type=l,
columns/T100rd/.style=column name=T100, column type=l,
columns/T50rd/.style=column name=T50, column type=l|,
columns/Evc/.style=column name=EAG, column type=l,
columns/T100vc/.style=column name=T100, column type=l,
columns/T50vc/.style=column name=T50, column type=l|,
columns/Emd/.style=column name=EAG, column type=l,
columns/T100md/.style=column name=T100, column type=l,
columns/T50md/.style=column name=T50, column type=l|,
columns/Ergas/.style=column name=EAG, column type=l,
columns/T100rgas/.style=column name=T100, column type=l,
columns/T50rgas/.style=column name=T50, column type=l|,
columns/Efdm/.style=column name=EAG, column type=l,
columns/T100fdm/.style=column name=T100, column type=l,
columns/T50fdm/.style=column name=T50, column type=l,

Physical parameters of galaxies in simulations (EAGLE: EAG; TNG100:T100; TNG50:T50). col:1 Bin number. col:2-4 Number of rotation curves in each bin. col:5-7 Disc radius. col:8-10 circular velocity computed at $R_{\mathrm{out}}$. col:11-13 stellar mass (includes bulge). col: 14-16 gas disc length. col:17-19 DM fraction. Here quantities are provided in $\log$ and their physical units are similar to the one given in Table <ref>.

Individual RC mass modelling: A few examples of full posterior distributions (MCMC output) of the estimated parameters ($\vec{\theta}$) for Burkert (red) and NFW (blue) halo. The vertical and horizontal gray dashed line on the posterior plots show the initial guess (or mean value in case of Gaussian prior) of parameter space $\vec{\theta}$.

Co-added RC mass modelling: A few examples of full posterior distributions (MCMC output) of the estimated parameters ($\vec{\theta}$) for Burkert (red) and NFW (blue) halo. The vertical and horizontal gray dashed line on the posterior plots show the initial guess (or mean value in case of Gaussian prior) of parameter space $\vec{\theta}$.