© ESO 2022

1. 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 (Öpik 1915; Kapteyn 1922; Jeans 1922; Oort 1932; Zwicky 1933; Smith 1936; Roberts & Rots 1973; Rubin et al. 1978; Bosma 1981; de Swart et al. 2017, and references therein). The astrophysical nature of such invisible matter can be inferred from its gravitational interaction on luminous matter (Ostriker & Peebles 1973; Persic & Salucci 1988; Persic et al. 1996; Navarro et al. 1997; Sofue & Rubin 2001; Salucci 2019, references therein) and on radiation which emerges at different cosmological scales (Sunyaev & Zeldovich 1969; Hu & Silk 1993; Boehm & Schaeffer 2005; Chluba & Sunyaev 2011; Dvorkin et al. 2014; Slatyer & Wu 2018; Ali-Haïmoud 2019, references therein). On the other hand, its particle nature, which is the center of the DM phenomenon, has not been characterized yet.

Recently, the DM phenomenon was uniquely framed into the dark energy driven cold dark matter (ΛCDM) scenario (Steigman & Turner 1985; Kolb & Turner 1990; Ellis 2018), where nonrelativistic DM can be described as a collisionless fluid, whose particle with a cross-section of about 10⁻²⁶ cm², only interacts gravitationally and very weakly with standard model particles. This scenario, which emerges very clearly cosmologically (Peebles 1993; Coles & Lucchin 2002), has been confirmed by large-scale observations (Einasto et al. 2017; Zaninetti 2018, and references therein) and can be investigated by appropriate simulations (Klypin 1996; Vogelsberger et al. 2020, and references therein) and semi-analytical modeling (Kauffmann et al. 1993; Somerville & Primack 1999; Lacey 2001; Guo et al. 2011; Somerville & Davé 2015, and references therein). As a consequence, a very well-defined perturbation power law spectrum and lack of dependencies on the initial conditions, yield 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, which is performed by various direct and indirect methods (e.g., Covi 2013; Gaskins 2016) and one of the major efforts of astro-particle physics, has been unsuccessful so far (Arcadi et al. 2018); (2) the main imprint of a collisionless particle is a DM halo with cuspy central DM density (ρ ∝ r⁻¹), which is not present in the majority of local spirals, dwarf disks, and low surface brightness galaxies (e.g. Palunas & Williams 2000; Salucci & Burkert 2000; de Blok et al. 2001; Karukes & Salucci 2017; Di Paolo et al. 2019). Moreover, local galaxies show that the physical properties (e.g., luminosity and disk scale length) of luminous matter are deeply connected with the structural properties of DM (e.g., its core radius and core density); readers can refer to Salucci (2019) and Li et al. (2020). This is very difficult to reconcile with collisionless 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, the DM halo evolution through cosmic time may also reflect the nature of a DM particle: in the ΛCDM scenario, the collisionless nature of a dark particle is itself imprinted in the dark halo evolution with time (Ellis 2018), which is well followed by hydrodynamical simulations (He 2020). Also in the other scenarios of DM such as warm DM, boson condensates, interacting DM, ultra light axions, etc., and with 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).

In the same vein, we must remind readers that, until recently, the determination of high-z DM halo density profiles was impossible: the kinematical measurements were lacking the sufficient sensitivity and spatial resolution. But, in the last decade, advanced use of integral field units (IFUs) in galaxy surveys has opened up possibilities for studying the spatially resolved kinematics of high-z galaxies, for example, surveys with the Multi-Unit Spectroscopic Explorer (MUSE: Bacon et al. 2010), K-band Multi-Object Spectrograph (KMOS: Sharples 2014), the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI: Eisenhauer et al. 2003), and Atacama Large Millimeter/submillimeter Array (ALMA: Cox 2016). A detailed review by Förster (2020, references therein) highlights the main findings of the aforementioned surveys. In the context of dynamical studies, Förster (2020) show that in using the resolved kinematics of galaxies at high-z, we can now obtain the rotation curves (RCs) of galaxies up to z ≈ 1 − 2.5. These RCs allow us to probe the baryon and DM content on galactic scales, as well as their distribution. In particular, Genzel et al. (2020) dynamically mass modeled the RCs of 41 star-forming galaxies (SFGs), and provided the first piece of evidence of a cored DM profile at high-z (in the majority of the sample). Moreover, in a more recent study, Bouché et al. (2022) dynamically mass-modeled the nine z ≈ 1 SFGs and they show that, similar to the z ∼ 1 sample of Genzel et al. (2020), a significant fraction of their sample also shows evidence of a cored DM profile. The study of Bouché et al. (2022) also shows that the DM “cored” feature is generally associated with galaxies with low stellar-to-halo mass, while “cusps” are found to be a function of strong star-formation activity. They conclude that this suggests feedback-induced core formation in the ΛCDM scenario. Now, in order to extend our understanding of high-z DM halos, in particular the cusp-core transformation, and to develop accurate models of galaxy evolution, we in this work focus on detailed mass modeling of the largest set of the currently available z ∼ 1 rotation curve of SFGs published and discussed in our previous works (Sharma et al. 2021a,b).

Briefly, in Sharma et al. (2021a), we used KMOS Redshift One Spectroscopic Survey (KROSS) data to derive the kinematics of a large sample of star-forming disk-like galaxies at z ∼ 1. We converted these observed kinematics into intrinsic rotation curves by correcting the issues related to high-z measurements, such as, beam-smearing¹ and pressure-support² conditions. In particular, we used 3D-forward modeling of datacubes using ^3DBarolo (Teodoro & Fraternali 2015), which provides somewhat better constraints on kinematics than 2D-kinematic modeling approaches. For pressure support, we employed the pressure gradient correction (PGC) method, which makes use of all of the information available in datacube and therefore provides better results than pressure support corrections that assume isotropic and constant velocity dispersion; for details we refer the interested reader to Sharma et al. (2021a, appendix). This allows us to have a large sample of more than 200 accurate RCs of disk-like galaxies at z ∼ 1 (i.e., 6.5 Gyr lookback time). This large amount of data can help us open a new portal to the nature of DM, as well as allow us to test the ΛCDM collisionless scenario.

In the follow-up study, Sharma et al. (2021b), we employed a halo model independent approach of DM determination using the RCs established in Sharma et al. (2021a), 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 (Tacconi et al. 2018; Lagos et al. 2011), and the dynamical mass (M_dyn = G⁻¹V²R) from the rotation curves, evaluated at a different scale radius (disk radius: 0.59R_e, optical radius: 1.89R_e, and outer radius: 2.95R_e ). These measurements provide us with DM fractions (f_DM = 1 − M_bar/M_dyn) at different radii. This work shows that the galaxies are DM dominated from 1 − 3R_e. The majority (> 72%) of star-forming galaxies in our sample have outer disks (∼5 − 10 kpc) dominated by DM, which agrees well with local SFGs.

In the present work, we investigate the kinematics of our same star-forming disk-like galaxies from z ∼ 1 and dynamically mass-model their RCs. The latter is a model-dependent approach for deriving the physical quantities such as the DM fraction, stellar masses, and effective radii. Thus, comparing the two studies, Sharma et al. (2021b) and the current work allowed us to check the robustness of our dynamical models. Through the mass-modeling techniques, we also determined the structural properties of DM halos and checked the scaling relations between dark and baryonic matter, which are observed in local star-forming disk galaxies. In the end, we compared our work with the recent study of Genzel et al. (2020), in which the authors dynamically mass model the galaxy kinematics, which was derived from the 1D software Slit Extraction Method. The latter method is known to bias the kinematic extraction (i.e., RCs), as shown in Teodoro & Fraternali (2015) and Di Teodoro et al. (2016). However, a complementary study of the work by Genzel et al. (2020) was done by Price et al. (2021) using 3D-forward kinematic models; however, they obtained the same results as Genzel et al. (2020). Besides this, both studies corrected the pressure support under the assumption of isotropic and constant velocity dispersion using Burkert et al. (2010) method, which underestimates the rotational velocity, as shown in Sharma et al. (2021a). Despite the significant discrepancies between our work and the aforementioned studies, we compared them in terms of the DM fraction and its structural properties.

This article is organized as follows: Sect. 2 contains a brief discussion on the data employed. In Sect. 3 we provide the mass models to disentangle the baryons and DM component from RCs and describe their fitting techniques. In Sect. 4 we show the results and critically analyze them. In Sect. 5 we compare our results with current state-of-the-art simulations. Finally in Sects. 6 and 7, we discuss the results and provide a summary with future prospectives. Throughout the analysis we assumed a flat ΛCDM cosmology with Ω_m, 0 = 0.27, Ω_Λ, 0 = 0.73, and H₀ = 70 km s⁻¹.

2. Data

We studied the KROSS parent sample previously studied by Stott et al. (2016) and others (Harrison et al. 2017; Tiley et al. 2019a,b; Johnson et al. 2018). This sample contains 586 SFGs from z ∼ 1. We selected 344 objects out of 586, on the basis of an integrated Hα flux cut (F_Hα > 2 × 10⁻¹⁷ [erg s⁻¹ cm⁻²]) and inclination angle (25° ≤θ_i ≤ 75°). The chosen flux and inclination cuts ensure sufficient signal-to-noise ratio (S/N) data and reduce the impact of extinction³ during the kinematic modeling. The intrinsic characteristic of the selected sample is the following (given with respect to Table A1 of Harrison et al. (2017, hereafter H17): (1)the AGN-flag is zero, that is there no evidence for an AGN contribution to the H_α emission-line profile; and (2)quality flag of objects is 1, 2, and 3, that is only H_α-emission line detected objects (S/N > 3) are selected. We adopted the values of requisite quantities, such as effective radii (R_e), stellar mass (M_*), and redshift (z) from Harrison et al. (2017) and Stott et al. (2016). Moreover, the reddening corrected Hα star-formation rate (SFR_Hα) of the sample were adopted from Sharma et al. (2021b). We remark that the parent sample was selected in such a way that it is biased against galaxies in a merging or interacting state, and thus should provide a representatives sample of main sequence SFGs at z ∼ 1 (see Stott et al. 2016).

As briefly mentioned above, the stellar masses and effective radii are two relevant physical properties of galaxies that we adopted from H17 and used throughout this work. In H17, stellar masses were calculated using a fixed mass-to-light ratio following M_* = Υ_H × 10^{−0.4 × (M_H − 4.71)}, where Υ_H and M_H are the mass-to-light ratio and absolute magnitude in the H band (rest frame), respectively. Here Υ_H = 0.2, which is the median value of our sample obtained using the HYPERZ (Bolzonella et al. 2000) spectral energy distribution (SED) fitting tool. These stellar masses are in good agreement with those inferred using a full SED-fitting approach (Tiley et al. 2019a), with no noticeable offset and only a 0.2 dex uncertainty (see Sharma et al. 2021b, appendix). The effective radii of the sample were measured from broadband images by deconvolving the point spread function (PSF) and semi-major axis of the aperture, which contains half of the total flux of a galaxy; for details, readers can refer to Harrison et al. (2017).

2.1. Individual RCs

To extract kinematics of a selected sample (344 SFGs), we exploited publicly available datacubes (data availability⁴), and the kinematics of these objects were re-derived from KROSS datacubes employing the 3D-BAROLO code (Teodoro & Fraternali 2015; Di Teodoro et al. 2016). The latter takes into account the beam smearing correction in 3D-space and provides moment maps, the stellar surface brightness profile, the RC, the dispersion curve (DC), along with the kinematic models. After running the 3D-Barolo, we assessed the quality of the objects, which is defined by taking the S/N into account, and limitation of 3D-Barolo (e.g., inclination). The best quality objects, that is to say quality-1 and 2, are the ones that have a S/N ≥ 3⁵ and an inclination 30° ≤θ_i ≤ 75°. Apart from this, our RCs were corrected for pressure gradients, which likely affect the kinematics of high-z galaxies. Out of the 344 KROSS objects, we chose to analyze only 256 rotation dominated objects, referred to as the Q12 sample, for details we refer the reader to Sharma et al. (2021a). In the following work, we use this Q12 sample to analyze the DM distribution.

2.2. Coadded RCs

In this work, we will investigate the individual RCs of the Q12 sample, as well as their coadded RCs (or stacked RCs). The process of coaddition, which we briefly recall here, is accurately described in Persic et al. (1996) and Salucci (2019). That is, the individual RCs were coadded 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 binned the galaxies radially per 2.5 kpc corresponding to the binning scale of the 3D-Barolo; for details, readers can refer to Sharma et al. (2021a). For the binning, we used the standard weighted mean statistic given by $\overline{X}=\frac{{\mathrm{\Sigma }}_{i=1}^n{x}_i\times {\mathit{w}}_i}{{\mathrm{\Sigma }}_{i=1}^n{\mathit{w}}_i}$, where, $x_i=\mathrm{data},{\mathit{w}}_i=1/{\mathrm{error}}^2,\mathrm{and}\overline{X}=\mathrm{binned}\mathrm{data}$. The errors on the binned data are the root mean square error, computed as ${\delta }_{\mathit{v}}^{r_i}=\sqrt{\sum ([{({\delta }_{{\mathit{v}}_i}^r)}^2+{({\sigma }_{\mathit{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}$ is scatter per radial bin. This process yields 16 coadded RCs, each one has a minimum of five independent data-points (see Fig. 1, and for more details on the binning process we refer the readers to Sharma et al. 2021a). While binning the RCs, we also coadded and binned the physical quantities (e.g., M_D, R_D) associated with individual objects per velocity bin using RMS statistics and the errors are standard errors. These binned parameters were used in RC modeling either as an 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 on top of them (e.g., ${\stackrel{\sim }{M}}_{\mathrm{D}}$) and coadded plus binned RCs are referred to as “coadded RCs”.

Fig. 1. Coadded and binned rotation curves. The individual and binned RCs are shown by blue thin curves and squares connected with a dark blue line, respectively. The central velocity (${\stackrel{\sim }{V}}_{\mathrm{c}}$) and the number of data points per bin are printed in the legend of each panel. For example, the notation ${\stackrel{\sim }{V}}_{\mathrm{c}}:82$ (15; 64) means that the central circular velocity of the bin is 82 km s−1 and that it contains 15 RCs and a total of 64 data-points. We note that here ${\stackrel{\sim }{V}}_{\mathrm{c}}$ is the circular velocity computed at ${\stackrel{\sim }{R}}_{\mathrm{out}}$ (=2.95 Re).

We remark that the main advantage of the coadded RC is that it allows for 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 high-z galaxies. Moreover, H_α emission is the only way to trace the kinematics of high-z galaxies in the optical (and near-infrared) wavelength range. At the same time we know that H_α originates from star-forming regions, that is the interstellar medium, which is turbulent at high-z; therefore, the individual RCs are generally noisy, that is to say artificially wiggled (Dyson & Williams 1980). In these mixed situations with poor data quality and limited observations, the coaddition of RCs acts as a rescue tool, providing a smooth circular velocity of the galaxies and reducing minor observational errors in the data. However, the coaddition of RCs also has some drawbacks, for example, the following can not be evident: (1) the diversity in the shape of the RCs of a given class of galaxies and (2) the information on noncircular motions due to a bar- and bulge-like component. The former and latter both are interlinked and arise due to a dominant triaxial halo which induces bi-symmetric flows in the gas component of a galaxy as reported in Marasco et al. (2018) and Oman et al. (2019). Despite these aspects, and taking the current situation of high-z observations into account where we have a lack of data and no resolution within 4 kpc, we propose to study coadded as well as individuals RCs.

Furthermore, before performing the coaddition of the RCs, we visually inspected each individual RC, and verified 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 a maximum dispersion of ∼0.1 dex, which is very similar to the study performed by Persic et al. (1996) for local SFGs. Thus, our method of coaddition is applicable and certainly necessary for high-z data, where we have moderate signal-to-noise, low resolution, and limited observations. An alternative way to avoid the abovementioned problems of high-z observations is “smooth RCs”, where we fit the observed RC with a parametric form, for example arctan (V_rot ∝ arctan(R)), then this best-fit RC can be used for mass modeling or coaddition. The main advantage of smoothed RCs is that while mass-modeling with Markov chain Monte Carlo (MCMC), they do not face the problem of nonconvergence. However, we have two important reasons not to use smoothed RCs. First, observed RCs are difficult to fit when they are wiggled (or declining), and thus they cannot be smoothed with simple parametric forms of RCs. Second, a statistical study of RCs, that is coaddition and binning performed in our work, leads to averaging RCs, which is similar to smoothing. Therefore, it is not advisable to smooth the RCs a priori, otherwise, twice the information may be lost. Moreover, our study focuses on the shape of the RCs (Sharma et al. 2021a), the DM fraction (Sharma et al. 2021b), the component separation (current work), and the derivation of the structural properties of the DM. These studies are very sensitive to the intrinsic shapes of the RCs. Therefore, we did not risk smoothing the RCs at outset but used the observed RCs.

3. Methodology

The 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 DM. Therefore, V_c(R) can be described as a sum in quadrature of the following:

$$\begin{array}{c}\hfill V_{\mathrm{c}}^2(R)=V_{\mathrm{D}}^2(R)+V_{\mathrm{bulge}}^2(R)+V_{\mathrm{H}2}^2(R)+V_{\mathrm{HI}}^2(R)+V_{\mathrm{DM}}^2(R),\end{array}$$(1)

where V_D(R),V_bulge(R),V_H2(R), and V_HI(R) are the circular velocity profiles of the stellar disk, stellar bulge, molecular (H2) and atomic (HI) gas disks, together providing the baryon’s contribution to the circular velocity and V_DM(R) gives the DM contribution. In the following sections, we present an approach to disentangle the aforementioned contributions from observed RCs. We note that in Sharma et al. (2021a), we corrected the observed RCs for inclination and pressure support, that is RCs are very close to intrinsic circular velocity curves.

3.1. Mass models

Here, we explain the various components of our mass modeling.

3.1.1. Stellar and gas disk

The stellar and gaseous components of local star-forming disks, dwarf-disks and low surface brightness galaxies are well described by the exponential disk model of Freeman (1970) (e.g., Persic et al. 1996; Karukes & Salucci 2017; Di Paolo et al. 2019). Here, we assume that the distribution of stars and gas in z ∼ 1 star-forming disk galaxies are also in exponential disks, so that the surface densities are as follows:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\_}(R)\propto \frac{M_{\_}}{R_{\mathrm{scale}}}exp(-R/R_{\mathrm{scale}}),\end{array}$$(2)

where M_{_} and R_scale are the total mass and the scale length of the different components (stars, H2, and HI), respectively. We note that stellar mass in the disk (without a bulge) is denoted by M_D, while the contribution of the total stellar mass is represented by M_star = M_D + M_bulge. Given the density distribution of the stars and the gas, their contribution to the circular velocity of the disk can be expressed as follows:

$$\begin{array}{c}\hfill V_{\_}^2(R)=\frac{1}{2}\left(\frac{G{M}_{\_}}{R_{\mathrm{scale}}}\right)(x^2)[I_0{K}_0-I_1{K}_1],\end{array}$$(3)

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 (cf. Persic et al. 1996; Karukes & Salucci 2017).

3.1.2. Bulge component

First, we would like to emphasize that our RCs are not well resolved within the 5 kpc where we expect the bulge to contribute significantly, that is, 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{array}{c}\hfill V_{\mathrm{bulge}}^2(R)=\frac{G{M}_{\mathrm{bulge}}}{R},\end{array}$$(4)

where, M_bulge is the total bulge mass, and it can be estimated using bulge-to-total ratio B/T = M_bulge/M_star.

3.1.3. DM halo component

In this work, we test two widely used DM halo models, namely the profiles of Burkert (Burkert 1995) and NFW (Navarro et al. 1996a).

Burkert halo. Local spirals and low surface brightness galaxies suggest existence of central DM cores, that is ρ_inner ∝ const (Salucci & Burkert 2000; de Blok et al. 2001). This observed DM-profile is well fitted by the Burkert (1995) halo, which possesses a double power law in the DM-density, that is, at small radii ρ ∝ R⁰ and at larger radii ρ ∝ R⁻³. Such a DM-density distribution can be expressed as follows:

$$\begin{array}{c}\hfill \rho (R)=\frac{{\rho }_0}{(1+\frac{R}{r_0})(1+\frac{R^2}{r_0^2}),}\end{array}$$(5)

where ρ₀ and r₀ are the central DM core density and core radius, respectively. Assuming spherical symmetry, the mass profile of the Burkert DM halo can be expressed as follows:

$$\begin{array}{c}\hfill M_{\mathrm{DM}}^{\mathrm{Burk}}(R)=4\pi {\rho }_0{r}_0^3[ln(1+\frac{R}{r_0})-arctan(\frac{R}{r_0})+0.5ln(1+\frac{R^2}{r_0^2})].\end{array}$$(6)

NFW halo. In the standard ΛCDM paradigm, the current cosmological simulations predict a cuspy DM distribution in the center, that is, ρ_inner ∝ R⁻¹. This type of DM profile is well approximated by the NFW halo, which is again a double power law; however, it has ρ ∝ R⁻¹ at small radii and ρ ∝ R⁻³ at larger radii. Such a DM-density distribution can be expressed as:

$$\begin{array}{c}\hfill \rho (R)=\frac{{\rho }_{\mathrm{s}}}{(\frac{R}{r_{\mathrm{s}}})(1+\frac{R}{r_{\mathrm{s}}}{)}^2},\end{array}$$(7)

where ρ_s and r_s are the characteristic density and scale radius of the DM distribution, respectively. Assuming spherical symmetry, the mass profile of the NFW DM halo is as follows:

$$\begin{array}{c}\hfill M_{\mathrm{DM}}^{\mathrm{NFW}}(R)=4\pi {\rho }_{\mathrm{s}}{r}_{\mathrm{s}}^3[ln(1+\frac{R}{r_{\mathrm{s}}})-\frac{\frac{R}{r_{\mathrm{s}}}}{1+\frac{R}{r_{\mathrm{s}}}}].\end{array}$$(8)

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

$$\begin{array}{c}\hfill V_{\mathrm{DM}}^2(R)=G\frac{M_{\mathrm{DM}}(<R)}{R},\end{array}$$(9)

where M_DM(< R) is the enclosed DM mass within radius R given in Eqs. (6) and (8). The above formalism allowed us to compute the various structural properties of the DM halo, namely r_0/s and ρ_0/s.⁶ For the purpose of this work, we also define the virial radius (R_vir), virial mass (M_vir), and concentration (c_vir) parameter of the halo. The virial radius is the distance within which the average halo density is 200 times the critical density (ρ_crit = 9.3 × 10⁻³⁰ gm cm⁻³) of the universe, it can be easily computed numerically. The virial mass is defined as follows:

$$\begin{array}{c}\hfill M_{\mathrm{vir}}=\frac{4}{3}\pi R_{\mathrm{vir}}^{3}200{\rho }_{\mathrm{crit}}.\end{array}$$(10)

In the case of the NFW halo, the ratio of the virial radius to scale radius gives the halo concentration, c_vir = R_vir/r_s. Furthermore, once the various components of the RC are disentangled, then, given the information on V_DM(R) and V_c(R), one can also estimate the following:

$$\begin{array}{c}\hfill f_{\mathrm{DM}}(<R)=\frac{V_{\mathrm{DM}}^2(R)}{V_{\mathrm{c}}^2(R)},\end{array}$$(11)

which gives us the DM fraction within the radius R.

3.2. Rotation curve modeling

To perform the mass modeling of RCs, we defined a χ² test statistic on our observed kinematics V_o(R) and modeled the kinematics V_m(R) as follows:

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{kin}}\equiv {\displaystyle \prod _{k=1}^{\mathrm{N}}\frac{1}{\sqrt{2\pi }\mathrm{\Delta }{V}_{\mathrm{o}}(R_k)}exp[-\frac{1}{2}(\frac{V_{\mathrm{o}}(R_k)-V_{\mathrm{m}}(R_k)}{\mathrm{\Delta }{V}_{\mathrm{o}}(R_k)}{)}^2].}\end{array}$$(12)

For each bin k ≤ N, we compared the theoretical prediction and the observed dataset by including the observational uncertainty ΔV_o(R_k) at each step. Our fitting procedure is based on Bayes’ theorem:

$$\begin{array}{c}\hfill \mathcal{P}(\mathrm{Model}(\mathit{\theta })|\mathrm{data})\propto {\mathcal{L}}_{\mathrm{kin}}(\mathrm{data}|\mathrm{Model}(\mathit{\theta })){\mathcal{P}}_0(\mathrm{Model}(\mathit{\theta })),\end{array}$$(13)

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 θ, with the likelihood function reported in Eq. (12), up to the overall normalization defining the so-called evidence of the model. The general model of the RC under the Bayesian inference is defined by five parameters θ = {M_D, R_D, M_bulge, r_0/s, ρ_0/s}. Here, three parameters, namely M_D, R_D, and M_bulge, constrain the baryons, and two parameters, namely r_0/s and ρ_0/s, describe the DM halo. The details of these parameters in various cases (studied in this work) of mass-modeling is given in Table 1. By employing Eqs. (12) and (13), we performed the MCMC analysis exploiting the RC dataset discussed in Sect. 2. To perform the MCMC, we made use of the emcee package (Foreman-Mackey et al. 2013). For each object, we allowed 100 walkers to evolve for 10 000 steps, which gave us 10⁶ samples for analysis. We removed the first half of the Markov chains to account for the burn-in period and the remaining were used to draw the posteriors of parameter space θ. Lastly, we tested our mass modeling technique on synthetic RCs. Those are modeled as in Sect. 3.1 using the realistic values for the parameter space θ. In Appendix B, we show that our MCMC approach provides an accurate fit to these synthetic RCs and that the retrieved parameter space θ is consistent, within the 1σ uncertainty, with the input parameters.

3.3. Fitting details

As we briefly mentioned above, in the fitting procedure, we have five free parameters, given by θ = {M_D, R_D, M_bulge, r_0/s, ρ_0/s}, where three parameters, namely M_D, R_D, and M_bulge, constrain the baryons and two parameters, namely r_0/s and ρ_0/s describe the DM in the case of the Burkert and NFW halo. The exploration range of each parameter is motivated by the mass modeling of local galaxies (see, Persic et al. 1996; Karukes & Salucci 2017; Lapi et al. 2018; Di Paolo et al. 2019), given in Table 1. As mentioned above, we performed RC modeling in two cases: case 1) individual RCs, and case 2) coadded RCs. In both cases, an initial guess on the parameters is the same and given as follows θ = {10, 0.3, 8.5, 1, −23.5} (in dex). The prior on the DM parameters were kept flat in both cases. On the other hand, the baryonic parameters have different conditions for prior in different cases, each motivated by the need to treat the data robustly and without biases. In the case 1, we used the Gaussian prior in the log scale for M_D, R_D, and M_bulge with a dispersion (σ) of 0.25 dex. In the log space, the unit is dex. Where the central values for M_D = 0.9M_star, which at first approximation excludes the bulge. The values of M_star and R_D are taken from Harrison et al. (2017). The bulge mass M_bulge is one-tenth of the stellar mass, which is motivated from local SFGs (for reference see Bruce et al. 2014; Morselli et al. 2017). In the case 2, we allowed the Gaussian prior only on ${\stackrel{\sim }{R}}_{\mathrm{D}}$, and ${\stackrel{\sim }{M}}_{\mathrm{D}}$, having a dispersion of 0.25 dex around the central value, and all other parameters kept a flat prior. In this case, central values are the quantities binned in the coaddition process of RCs, as discussed in Sect. 2.2. Similar to the case 1, central value of ${\stackrel{\sim }{M}}_{\mathrm{D}}$ is $0.9{\stackrel{\sim }{M}}_{\mathrm{star}}$.

We note that our modeling approach also requires other parameters in the fitting process, namely M_H2, M_HI, R_H2, and R_HI (see Eq. (1)). In particular, M_H2 and M_HI were calculated in each MCMC run for a given star formation rate and redshift (M_H2, HI(SFR, M_D, z)) using the scaling relations given by Tacconi et al. (2018) and Lagos et al. (2011), respectively. The molecular gas scale length R_H2 was determined from the Hα surface brightness distribution, since ionized gas around a star is surrounded by molecular gas clouds. For a detailed calculation of M_H2, M_HI and R_H2, we refer the reader to Sharma et al. (2021b). 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 all Fig. 5) see also (Leroy et al. 2008; Cormier et al. 2016). Therefore, we assume R_HI = 2 × 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. Furthermore, it is important to note that these additional parameters, such as the star-formation rate (SFR), M_gas, and R_gas, were fixed in the modeling procedure to the best of our knowledge.

While performing our mass modeling on individual RCs (i.e., Case-1), we visually (qualitatively) inspected the individual RCs and each of their individual fittings. We found that there are 28 perturbed RCs that we could not model, that is, we mass-modeled 228 individual RCs⁷. After the mass modeling 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_D). This tells us that the fitting of these 92 objects is prior-driven, so we discarded them from the individual analysis. Thus, we ended up analysing only 136 individual RCs. It is important to note, however, that in case 2, we used up all of the 228 RCs, making 16 coadded RCs (see Sect. 2.2), which were mass modeled using the aforementioned prescription.

4. Results and analysis

As described in Sects. 3.1 and 3.2, by inputting the individual RCs (case 1) and the coadded RCs (case 2), we constrained the parameter space θ and disentangled the various components of each RC. With this information, we estimated the quantities R_vir, M_vir, c_vir, and f_DM. In the following section, we present these results.

4.1. Mass modeled RCs

Here, we present the results of mass modeling in case 1 and case 2. Some examples of their MCMC outputs (posteriors) and best fit to the RCs are shown in Figs. 2 and 3, respectively. We note that in both cases, the Burkert and NFW halo models provide similarly good fits (in terms of minimum ${\chi }_{\text{red}}^2$) to the data. The rationale 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 just by looking at the rotation curve fitting.

Fig. 2. Individual RC mass modeling. Row 1: posterior distributions (MCMC output) of the following estimated parameters: r0 and ρ0 for Burkert (red) and rs and ρs for the NFW (blue) halo. The vertical and horizontal gray dashed lines show the initial guess for parameters r0/s and ρ0/s. The full posterior distribution of parameter space θ of these example cases are shown in Fig. D.1. Rows 2 and 3: best fit to the RCs and their velocity decomposition for Burkert and NFW halos, respectively. In the 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σ 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 the gray shaded area represents the unresolved region of the RC. The best-fit value of parameter space θ is given in the legend of the figures. All the parameters are printed in log, having physical units as follows: MD [M⊙], RD [kpc], Mbulge [M⊙], r0/s [kpc], and ρ0/s [gm cm−3].

Fig. 3. Coadded RC mass modeling. Row 1: posterior distributions (MCMC output) of the following estimated parameters: ${\stackrel{\sim }{r}}_0$ and ${\stackrel{\sim }{\rho }}_0$ for Burkert (red) and ${\stackrel{\sim }{r}}_s$ and ${\stackrel{\sim }{\rho }}_s$ for the NFW (blue) halo. The vertical and horizontal gray dashed lines show the initial guess for the parameters ${\stackrel{\sim }{r}}_{0/s}$, ${\stackrel{\sim }{\rho }}_{0/s}$, and ${\stackrel{\sim }{M}}_{\mathrm{bulge}}$. The full posterior distribution of parameter space θ of these example cases are shown in Fig. D.2. Rows 2 and 3: best fit to the RCs and their velocity decomposition for the Burkert and NFW halos, respectively. In the 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σ 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 the gray shaded area represents the unresolved region of the RC. The best-fit value of parameter space θ is given in the legend of the figures. All the parameters are printed in log, having physical units as follows: ${\stackrel{\sim }{M}}_{\mathrm{D}}[M_{\odot }]$, ${\stackrel{\sim }{R}}_{\mathrm{D}}[\mathrm{kpc}]$, ${\stackrel{\sim }{M}}_{\mathrm{bulge}}[{\mathrm{M}}_{\odot }]$, ${\stackrel{\sim }{r}}_{0/s}[\mathrm{kpc}]$, and ${\stackrel{\sim }{\rho }}_{0/s}[\mathrm{gm}{\mathrm{cm}}^{-3}]$.

Next, we focused on comparing the physical properties of galaxies, such as the total stellar and bulge masses. In particular, we compared the physical properties derived from case 2 because this case gives us a smooth distribution of RCs that reduce the effect of minor observational errors (perturbations caused by bad data points), which make the mass modeling more impeccable than in case 1. In the left panel of Fig. 4, we compare the dynamically derived stellar masses from the 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 dex 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.

Fig. 4. Dynamically derived integrated quantities. Left panel: comparison of the photometric and dynamical stellar masses in the Burkert (blue) and NFW (pink) halo profiles for coadded RCs. The black solid line shows the one-to-one relations. We note that the dynamical stellar masses obtained assuming 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 dex. Right panel: dynamically derived bulge masses versus stellar masses, in the case of Burkert (blue) and NFW (pink) halo profiles for coadded RCs. We compared our results of the same redshift range (z = 0.6 − 1.2) and a similar stellar mass sample with a pioneering study by Genzel et al. (2020, red filled circles), and we found that in our sample B/T is never higher than ∼0.3, which is different from what is reported by Genzel et al. (2020).

Moreover, in both cases (with both halo models), we found it difficult to constrain the bulge mass, and often ended up with a broad distribution. However, the bulge is never heavier than 0.3 × M_star (see right panel of Fig. 4) and, therefore, it does not significantly affect the inferred DM fraction (see Fig. 5). In fact, the DM has emerged as a dominant quantity, especially in the outskirts⁸, which is very similar to the local SFGs (Persic et al. 1996; Sharma et al. 2021b). We also found that for low mass objects (M_D < 10^9.5 M_⊙), the atomic gas kinematics (HI) dominate the baryonic component; whereas for massive objects, the stellar kinematics prevail.

Fig. 5. DM fraction within Rout for the Burkert and NFW halos, on the left and right, respectively. The color coding in both panels are the same for the individual RCs (blue filled circles), coadded RCs (blue hexagons), local SFGs (coral red circles), and high-z SFGs (coral red hexagons). The size of the makers in local SFGs includes their errors. The blue shaded area represents the model-independent DM fraction study that we previously conducted in Sharma et al. (2021b).

4.2. DM fraction

In Fig. 5 we show the results of the DM fraction within R_out (=2.95R_e) inferred from our mass modeled RCs, in case 1 and case 2, for the Burkert and NFW halo. Firstly, we noticed that both halo models in both cases have shown that the majority (∼70%) of galaxies are DM-dominated in the outskirts. Secondly, our current results are in very good agreement with our previous study (Sharma et al. 2021b), where we determined the DM fraction using a halo model-independent approach. This confirms that the majority of SFGs at z ∼ 1 have DM-dominated outer disks (up to 3 × R_e).

In the case of individual RCs, in both halo models, nearly 20% of objects appear to be deficient in DM (i.e., f_DM(< R_out) < 50%). We also note that the DM fraction derived from the individual RC modeling has a higher uncertainty in both halo models, while it is significantly reduced in the case of the coadded RCs. This could be due to the low resolution in the data (nonsmooth V(R)), which is penalized by the large uncertainty. On the other hand, coadding RCs statistically reduces the effect of the random noise introduced by the bad data points and provides us with a smooth RC, which can be easily dynamically modeled, and hence yields a relatively small uncertainty on the measured parameters.

As we discussed in Sect. 4.1, dynamical modeling assuming a NFW halo profile tends to suppress stellar mass (see left panel of Fig. 4). Therefore, the DM fraction obtained in this case is higher in comparison to the Burkert halo case (see right panel of Fig. 5). The latter argument is clearly visible from a coadded RC (blue hexagons) analysis. This suggests that dynamical modeling with the NFW halo requires precise baryon information, while the Burkert halo does not and it is still able to constrain the physical properties of galaxies, for example stellar masses. This is the reason that the f_DM − M_* relation obtained with the Burkert halo matches that inferred in Sharma et al. (2021b) very well for which a halo model-independent approach was used (blue shaded region in Fig. 5). However, this cannot be a strong piece of evidence to discard one model over the other because measurements from NFW models are still in good agreement within the 1σ uncertainty.

4.3. Structural properties of DM

With the information on the DM halo parameters r_0/s and ρ_0/s, we derived the physical properties of the galaxies, such as R_vir and M_vir, and the concentration in the case of the NFW halo. Using the latter quantities, we compared the c_vir − M_vir relation of our sample with DM-only simulation (Bullock et al. 2001) and the NIHAO hydrodynamical simulation (Wang et al. 2015), as shown in Fig. 6. We found that the concentration in z ∼ 1 SFGs goes in the opposite direction, which challenges the simulations. In other words, both the DM-only and hydro simulations predict the halo concentration which decreases with increasing redshift (see red and blue dashed lines in Fig. 6). The observations show a relatively high concentration at z ∼ 1, which is somewhat in tension with current simulations, although they remain consistent with them within the 1σ uncertainty. We remark that the concentration was derived as c_vir = R_vir/r_s, where the determination of R_vir requires a huge extrapolation. Therefore, the concentration in this work is not as accurate as it can be derived from HI RCs, where data are available over more than a 50 kpc radius. Thus, although the NFW halo gives conflicting results with respect to the simulations, we cannot falsify its presence in nature.

Fig. 6. Concentration vs. virial mass relation for individual and coadded RCs (or galaxies) in blue circles and hexagons, respectively. The blue and red, solid and dashed lines represent the c − Mvir relation of DM-only simulations (Bullock et al. 2001) and the NIHAO suite of hydrodynamical cosmological simulations (Wang et al. 2015; Freundlich et al. 2020), for z = 0 and z ∼ 1, respectively.

On the other hand, the relation between the core density and core radius, associated with the Burkert profile, emerge very similar to the local SFGs, with a similar slope, as shown in Fig. 7. However, we found 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 ∼ 1 are smaller by a factor of 0.3 dex and denser by 1.5 dex.

Fig. 7. DM core density vs. core radius relation of local and high-z galaxies. The blue circles and hexagons represent the z ∼ 1 individual and coadded RCs, respectively. The coral red circles and stars represent the local star-forming and low surface brightness galaxies (Persic et al. 1996; Di Paolo et al. 2019), respectively. The brown dashed line is the best fit ${\rho }_0\propto r_0^{\alpha }$ to z ∼ 1 data with a slope α = −1.54, and the orange shaded area represents the 1 − 3σ intrinsic scatter in the relation.

5. Comparison with simulations

In this section, we compare the mass modeling 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 Appendix A. In particular, we compare our observed coadded RCs with the simulated coadded circular velocity curves. Briefly, in simulations we know the exact total circular velocity of the galaxy as a function of radius (V_c(R)), and its various components, such as V_star(R), V_DM(R), and V_gas(R). From observations we only know the V_c(R), which we mass modeled to disentangle the various components (see Sect. 3.1). That is to say, in both cases, observations and simulations, we know the V_c(R) and V_DM(R). This allowed us to compute the DM fraction of galaxies as a function of radius using Eq. (11). In Fig. 8 we show the results of f_DM(< R) in the aforementioned simulations. In particular, we show 16 velocity bins, in which the first four bins have a few galaxies (in particular EAGLE simulations), and the remaining 12 bins have enough simulated data (in all simulations) to bin and obtain the coadded RCs; for details, readers can refer to Table C.2. It is noticeable that all simulations show DM-dominated galaxies from 5 kpc outward. 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 ∼ 1, reported previously in Sharma et al. (2021b). We remind the reader that galaxies within effective radii are known to have a lower DM fraction because baryon dominates the inner region of galaxies (Persic et al. 1996, and references therein). If one has to look for DM, it should be found in the outskirts (> 2 × R_e), as clearly shown in Fig. 8.

Fig. 8. DM fraction as a function of the radius in various velocity bins of observed and simulated galaxies. The central value of fDM(< 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σ error in each measurement is represented by a shaded area; whereas in the observations, it is a dashed line following the same color code. The gray shaded area in each panel represents the unresolved region in the observations. The horizontal red dotted-dashed line in each panel shows the 51% DM fraction. The central value of the circular velocity of each bin is printed at the top of each panel, where ${\stackrel{\sim }{V}}_c^{\mathrm{SIM}}$ and ${\stackrel{\sim }{V}}_c^{\mathrm{OBS}}$ specify the simulations and observations, respectively.

To compare the simulations with the observations, we superimposed the f_DM(< R) of the observed coadded RCs in Fig. 8, 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. It is important to notice that we compare the observed and simulated coadded RCs only when $V_{\mathrm{c}}^{\mathrm{SIM}}\approx V_{\mathrm{c}}^{\mathrm{OBS}}$, as shown at the top of each panel. First, we noticed that in observations, we have roughly three coadded RCs, corresponding to log(M_star [M_⊙]) ≈ 9.2, 9.4, 10.0, with a low DM fraction (i.e., f_DM < 50%) in the outskirts (for Burkert: $V_{\mathrm{c}}^{\mathrm{Obs}}\sim 83,114{\mathrm{km}\mathrm{s}}^{-1}$ and NFW $V_{\mathrm{c}}^{\mathrm{Obs}}\sim 83,114,127{\mathrm{km}\mathrm{s}}^{-1}$). We found that these RCs have a gas-dominated outer disk, which we could not find in the simulations. This could be either due to limitations of simulations, or scaling relations that we used in gas mass estimates. Despite very low velocity bins (V_c < 130 km s⁻¹) for both Burkert and NFW halos, we noticed DM-dominated outer disks. Second, the NFW and Burkert halo results at R > 10 kpc are very similar in intermediate velocity bins ($V_{\mathrm{c}}^{\mathrm{SIM}}\approx 130-170\mathrm{km}{\mathrm{s}}^{-1}$) and they agree well with the simulations within a 1σ uncertainty. However, when we go toward the higher velocity range, the NFW halo starts to show a deviation from the simulations. In particular, the Burkert halo and the simulations have almost the same DM fraction between the 10–15 kpc radius, which is different from the NFW halo. In short, the Burkert profile seems to represent the simulations more accurately.

To observe the largest differences, we emphasize that a reader should look at the DM fraction between the 5–10 kpc radius in Fig. 8, where we can clearly notice that the NFW halo allows for more DM content than for the Burkert halo, which we do not even see in the simulations. In other words, the DM fraction derived from the Burkert halo model is close to the one seen in hydrodynamical simulations for the majority of circular velocity bins, while NFW yields a relatively high DM fraction. The latter could be a consequence of suppressing the stellar mass which we have already noticed from the mass modeling of observations in Sects. 4.1 and 4.2, and Figs. 4 and 5.

6. Discussion

To fit the individual RCs through a mass modeling approach, one needs to use the informative prior on the baryons (M_D, M_bulge, and R_D), otherwise the parameter space θ is highly degenerated. In the case of coadded RCs, on the other hand, we could achieve a good fit only by specifying two baryonic parameters, which are the stellar disk scale length and stellar mass in this work. This case also allowed us to adopt the flat prior on bulge mass, which is an unknown quantity (and unresolved), and thus we can determine it dynamically (see Sect. 4.1 and Table 1). Based on RC fitting and structural properties of DM, we cannot robustly validate or invalidate the NFW halo model, but tension exist. For instance, a galaxy cannot have a stellar mass lower than what is measured from photometry, and if so, it misleads the mass budget of a galaxy. The latter happens in the case of the NFW halo model fitting, where mass modeling yields stellar masses that are lower than that inferred from the observed photometry. This results in a high DM fraction that can be seen in the right panel of Fig. 5, which we do not see when using the Burkert halo, nor in the simulations (EAGLE, TNG100, and TNG50; see Fig. 8). Providing the various results and comparisons of the NFW and Burkert halo in the case of individual and coadded RCs, we found that our current sample favors the Burkert halo model and also agrees quite well with the simulations. Moreover, the recent studies by Genzel et al. (2020) and Bouché et al. (2022) also found a cored profile in majority of their SFGs at z ∼ 1. Therefore, we focus on discussing the Burkert halo results in what follows.

We begin by referring to our earlier study of the DM fraction (Sharma et al. 2021b), in which we used a halo model-independent approach to determine the f_DM from RCs, and showed that the galaxies are DM-dominated until the last point of observation (1 − 3 R_e). We found that our earlier finding agrees very well with our current results of the DM fraction obtained via mass modeling (see the left panel of Fig. 5). These results are also consistent with the previously studied local and high-z SFGs, indicated by coral red circles and hexagons, respectively (see Fig. 5). In this work, we do not compare the DM fraction within R_e because it is beyond the scope of this paper. However, it is done in our previous work (Sharma et al. 2021b, see Fig. 5), where we have compared f_DM(< R_e) of our sample with the findings of Genzel et al. (2020). We have shown that ∼25% of their z ∼ 1 sample agrees well with our work, and the rest of the sample shows a low DM fraction (i.e., baryon-dominated disk) as their sample is biased toward massive galaxies. Now, with dynamical mass-modeling we refined our previous work, which allowed us to compute the DM fraction as a function of the radius of a galaxy (see Fig. 8). As one can notice, the majority of the sample shows DM-dominated outer disks (from 5–15 kpc, shown by a horizontal orange dotted-dashed line), which falsifies the general remark of a low and modest DM fraction in z ∼ 1 SFGs, made in Genzel et al. (2020). Meanwhile, our new results of f_DM(< R) match with the simulations (see Fig. 8). Furthermore, we would like to mention that the lower velocity end (V_c < 130 km s⁻¹) of our pilot study is less understood, firstly because we have very few (one to three) galaxies in simulations, which do not allow us to perform an accurate statistical analysis and, secondly, scaling relations that allow us to estimate the gas content from observed galaxies are likely not applicable on small systems (i.e., V_c < 130 km s⁻¹). Therefore, in the future, we will improve upon these two aspects using ALMA observations, upcoming HI-surveys, and by enlarging the simulated sample.

In one of our previous works (Sharma et al. 2021a), we have shown that the progenitors of present-day SFGs at z ∼ 1 have similar rotation curves. Given the similarity in the RCs of z ∼ 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 ∼ 1 to z = 0, while the stellar mass distribution (traced by the light profile) evolves. To bring together our previous and current findings, we compared the structural properties of the DM halo at z ∼ 1 and z ≈ 0 in the circular velocity plane (Fig. 9). We would like to remark that the r₀ − V_c and ρ₀ − V_c relations of local galaxies are well fitted by second order polynomials (see upper panel Fig. 9), which are given by the following:

$$\begin{array}{c}\hfill log{r}_0^{z\sim 0}=1.63{(log{V}_{\mathrm{c}})}^2-5.24log{V}_{\mathrm{c}}+4.76,\end{array}$$(14)

$$\begin{array}{c}\hfill log{\rho }_0^{z\sim 0}=-1.99{(log{V}_{\mathrm{c}})}^2+6.67log{V}_{\mathrm{c}}+2.04.\end{array}$$(15)

Fig. 9. Structural properties of the DM halo at z ∼ 1. Top: DM core radius and core density as a function of the circular velocity shown in the upper and lower panels, respectively. The individual and coadded RC measurements are represented by small blue circles and hexagons, respectively. The comparison sample of local LSBs (Di Paolo et al. 2019) and SFGs (Persic et al. 1996; Lapi et al. 2018) 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 r0 − Vc plane and below in the ρ0 − Vc plane with respect to the z ∼ 1 SFGs. This is due to the fact that DM cores of z ∼ 1 SFGs are smaller and denser than their local counterparts, indicating an evolution in the structure of DM halos over cosmic time. Bottom: ratio of DM core radius and density at z ≈ 0 to z ∼ 1 as a function of circular velocity. This plot shows the expansion factor of DM cores at z ∼ 1 with respect to local SFGs. The color codes for data points are similar to those in the top 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 a solid line shows the 1σ intrinsic scatter in the relations.

Therefore, to visualize the differences in the structural properties of DM at different epochs, we used Eqs. (14) and (15) as a benchmark and plotted the ratios: $\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 Fig. 9. This immediately shows us how the structural properties of DM differ at z ∼ 1 with respect to z ≈ 0. First, we find that the DM cores at z ∼ 1 are smaller than their local counterparts by a median factor of 0.3 dex. Second, the density of DM cores at z ∼ 1 is 1.5 dex higher than at z ≈ 0. It is also very clear that these two ratios correlate very well with the circular velocity of the galaxies. We refer to these ratios as the expansion factors of the core radius and density, which are both a function of the circular velocity and well fitted by the following:

$$\begin{array}{c}\hfill log(\frac{r_0^{z\sim 0}}{r_0^{z\sim 1}})=2.1{(log{V}_{\mathrm{c}})}^2-7.3log{V}_{\mathrm{c}}+6.6\lfloor \pm 0.6\mathrm{dex}\rceil ,\end{array}$$(16)

$$\begin{array}{c}\hfill log(\frac{{\rho }_0^{z\sim 0}}{{\rho }_0^{z\sim 1}})=-2.7{(log{V}_{\mathrm{c}})}^2+8.2log{V}_{\mathrm{c}}-5.8\lfloor \pm 0.9\mathrm{dex}\rceil .\end{array}$$(17)

The Eqs. (16) and (17) are accurate in relating the structural properties of DM to baryonic matter. Plotting these relations in the lower panel of Fig. 9 clearly shows us that at z ∼ 1, objects that have a higher circular velocity are denser in nature and more compact in size (in terms of DM cores), relative to z ≈ 0 SFGs. In other words, objects with a higher circular velocity expand their DM cores more over cosmic time, compared to objects with a lower circular velocity. This means that the circular velocity plays a crucial role in the galaxy evolution. We report that the intrinsic scatter in relation (16) and (17) is about 0.6 and 0.9 dex, respectively. In addition to 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 Genzel et al. (2020). This discrepancy is likely due to the different methods and corrections applied in deriving the intrinsic RCs, which certainly affect the determination of the halo parameters. For example, most of the Genzel et al. (2020) RCs are declining, that is they are baryon-dominated systems. Therefore, the mass-modeling of these RCs yields a larger core radius of the DM halo and consequently a low concentration parameter.

The discussion above 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 if theory is correct and that early DM halos are cuspy (for instance at z > 2), they must have changed over cosmic time under the influence of galactic processes (e.g., AGN and supernova feedback as well as dynamical friction), which has been previously suggested by many authors such as Navarro et al. (1996b), Pontzen & Governato (2012), Nipoti & Binney (2014), Read et al. (2016, 2018), and Lazar et al. (2020). Although, the difference in the size and density of the DM core between z ∼ 1 and z ≈ 0 is small, it is nevertheless remarkable and could have a very strong impact on the identification of the nature of DM and building the models for galaxy evolution. In short, our results suggest that the structure of DM halos expand over cosmic time. Based on these results, one can easily interpret that if the DM cores are denser at z ∼ 1, then the virialization scale (R_vir) 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 (Baker et al. 2018).

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. To confirm our findings, we still require high resolution and good S/N data in the inner region of the galaxies. That can possibly be obtained with the upcoming James Webb Space Telescope (JWST) and in future with the Extremely Large Telescope (ELT).

7. Summary and conclusions

In this work, we have studied the structural properties of the DM in z ∼ 1 star-forming disk-like galaxies. We first mass-modeled the individual and coadded rotation curves, assuming that the baryons are distributed in the Freeman (1970) disk and that DM is modeled with the Burkert and NFW halo profiles. Then the mass modeling results in the Burkert and NFW cases were critically examined, interpreted, and compared with current state-of-the-art galaxy simulations. We found that the results obtained using the Burkert halo fit well with the observations and also match with the simulations (see Figs. 2–5, and 8). Therefore, in the rest of our paper, we solely discussed measurements obtained with the Burkert halo profile. We found the following:

• Our mass modeling approach has been proven to be robust to disentangle the various components of RCs, and it was able to determine the 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 DM (for details see Sect. 4 and Appendix B).

• The DM fraction derived from the mass modeling outcomes (using Eq. (11)) agrees very well with the results of the halo-model independent DM fraction (obtained in Sharma et al. 2021b). Moreover, it is also similar to those predicted by different simulations: EAGLE, TNG100, and TNG50.

• The results of the DM structural properties (r₀ and ρ₀) show that the DM cores at z ∼ 1 are, on average, a factor of 0.3 dex smaller and 1.5 dex denser than their local counterparts. This suggests an expansion in the structure of DM halos over the last 6.5 Gyrs.

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

We have shown in our previous work that the progenitors of local star-forming galaxies at z ∼ 1 are similar in terms of total mass ($V_{\text{c}}^2(R)\propto M_{\text{tot}}$); they are DM-dominated objects, and their stellar mass distribution evolves with time as a function of their circular velocity (see, Sharma et al. 2021a,b). Now, this pilot study of dynamical mass modeling of high-z RCs shows that DM 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 DM, and consequently DM distribution changes. In other words, DM responds to baryonic processes (e.g., AGN and supernova feedback) over the timescale of a few gigayears. 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 and reviewed by Pontzen & Governato (2014, and references therein). This piece of evidence will have a strong implication for cross-checking simulations, modeling galaxy evolution, as well as determining the astrophysical and particle nature of DM.

For physical properties of galaxies, such as stellar mass, gas mass, SFR, velocity, we refer readers to the catalogs released with Sharma et al. (2021a) (Flat RCs catalog, Sharma 2021a) and Sharma et al. (2021b) (DM fraction catalog, Sharma 2021b). The measurements from coadded RC mass modeling are available in form of tables in Appendix C.

Acknowledgments

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 as well as discussion and 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).

References

Appendix A: Simulations

In addition to observations, we also analyzed recent state-of-the-art galaxy simulations, EAGLE (Schaye et al. 2015) and IllustrisTNG (Pillepich et al. 2018). For the purposes of this paper, we selected central galaxies from the simulations using a friends-of-friends halo finder approach (Davis et al. 1985) and the SUBFIND algorithm (Springel et al. 2001; Dolag et al. 2009). Gas and star particles are associated with the galaxy of their nearest neighboring DM particle if they belong to one, and if they are gravitationally bound to the object. The center of a galaxy is defined as the location of its particle and 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 3.1.3.

To match the observed and simulated galaxies, we further filtered the central galaxies using the following selection criteria: z = 0.85, M_star = 10⁸ − 10^11.6 M_⊙, SFR = 0.1 − 100 M_⊙yr⁻¹, R_e = 0.03 − 30 kpc, V_c = 75 − 300 km s⁻¹, and V_rot/σ > 1. Here, V_c is the maximum circular velocity of the galaxies within the radius of 100 kpc. The V_rot/σ is the ratio between the 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 A.1, where we show the main sequence of SFGs, the mass-size relation, and the velocity distribution. 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 (Whitaker et al. 2014; Donnari et al. 2019). 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 ∼ 1 in the ΛCDM universe.

Fig. A.1. Comparison of the physical parameters of the 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, the so-called main sequence of SFGs. Middle panel: Effective radii of galaxies as a function of stellar mass, the so-called mass-size relation. Right panel: Distribution of the maximum circular velocity of the simulated sample compared to Vout = Vc of the observations.

To compare the simulations with observations, we extracted the mass profile of simulated galaxies expressed as the circular velocity curve ($V_c(R)=\sqrt{GM/R}$) and its constituent components (i.e., stars, gas, and DM). We coadded the simulated circular velocity curves using a similar approach as used for observations, and binning was performed using the RMS statistic. Briefly, we divided our simulated galaxies into 20 velocity bins between 100-300 km s⁻¹ with a step size of 10 km s⁻¹. Then each bin was coadded and binned to obtain a single circular velocity curve. We would like to note that in the simulations, we know the exact total circular velocity of the galaxy, as well as its different components, namely stars (V_star), gas (V_H2 and V_HI), and DM (V_DM). Therefore, in the binning process, we binned the total circular velocity, as well as its various components. Moreover, we also binned the physical quantities associated with each galaxy, such as R_D and M_star. In Figure A.2, we show an example case of coadding and binning, as well as the final coadded circular velocity curves from all three simulations, and in Table C.2, we tabulate all of their relevant physical quantities. These circular velocity curve data were then used in the investigation of the DM fraction (see Section 4). We remark that the TNG100 and TNG50 simulations are the same and only differ in terms of resolution. However, we still notice a small dissimilarity 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 it is beyond the scope of this work and will be covered in our future work for which we plan to make the synthetic datacubes of simulated galaxies and determine the RCs using the 3D Barolo, as suggested by Oman (2017) and Oman et al. (2019).

Fig. A.2. Individual 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 the velocity bin is shown in the title of each plot.

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

Fig. A.4. Binned RCs of all velocity bins of galaxy simulations from first to last EAGLE, TNG100, and TNG50. The color maps give an idea of their increasing circular velocities.

Appendix B: Testing mass modeling approach

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

B.1. Synthetic RCs

To generate the synthetic RCs, we took the reverse approach, that is, we used the best-fitting parameters obtained from the mass modeling of the individual RCs, namely θ = [M_D, R_D, M_bulge, r_0, s, ρ_0, s], gaseous disk parameters = [M_H2, M_HI, R_H2, R_HI], and SFR. The redshift of the synthetic RCs was fixed to z = 0.85. Then, the synthetic RCs were generated using the baryon and DM models defined in Section 3.1, and in the end we assigned to each RC the uniform error of 5-25 km/sec typically observed in the true RCs. In this way, we obtained 228 individual synthetic RCs, for which we know the exact baryon and DM content. An example of such individual RCs is shown in the bottom panel of Figure B.1. We note that we generated synthetic RCs for both the Burkert and NFW halo profiles. We used these synthetic RCs to test our mass modeling approach.

Fig. B.1. Example of synthetic coadded RCs. Top panel: Observed coadded and binned RC, which is shown as a reference. Bottom panel: Coadded and binned synthetic RCs in the case of the Burkert and the NFW dark matter halo, which are on the left and right, respectively. The color code in all of the panels are the same and given as follows: the brown square represents the individual RCs and the blue stars represent the coadded RC.

B.2. Synthetic coadded RCs

Coaddition of synthetic RCs was performed using the same strategy applied on observed RCs, defined in Section 2.2. That is, each synthetic coadded RC contain 15 synthetic individual RCs; they were binned using weighted mean statistics. During the binning process, we also binned the parameters space θ and gas disk parameters, yielding ${\stackrel{\sim }{M}}_D,{\stackrel{\sim }{R}}_D,{\stackrel{\sim }{M}}_{\mathrm{bulge}},{\stackrel{\sim }{r}}_0,{\stackrel{\sim }{\rho }}_0,{\stackrel{\sim }{M}}_{H2},{\stackrel{\sim }{M}}_{\mathit{HI}},{\stackrel{\sim }{R}}_{H2},$ and ${\stackrel{\sim }{R}}_{\mathit{HI}}$ for synthetic coadded RCs, which we used as a central values in mass modeling (if required). An example of a synthetic coadded RC is shown in Figure B.1.

B.3. Mass modeling of synthetic coadded RCs

For the mass modeling of synthetic coadded RCs, we used the same technique as applied on observed coadded RCs, and described in Section 3.2. We examined the parameter space (θ) using the MCMC sampler. In the fitting process, we kept the Gaussian prior in the log scale for R_D and M_D with a 0.25 dex dispersion around their central values. For the rest of the parameters (θ), we used the flat prior in the log scale with ranges given in Table 1. The gas mass (atomic and molecular) depends on the scaling relations (Tacconi et al. 2018; Lagos et al. 2011) and it was 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 scale lengths of the gaseous disks (R_H2 and R_HI) are fixed.

Results: In the lower panel of Figure B.1, we show an example of synthetic individual and coadded RCs for bin number five (bin_5) in the case of Burkert and NFW halo profiles, which are on the left and right, respectively. These synthetic RCs are shown in comparison with ther observed coadded RCs (upper panel) of the same velocity bin. The binned values (i.e., the true values) corresponding to each coadded RC are given in the title of each RC plot. We allowed mass modeling of these RCs using the aforementioned procedure. In Figure B.2, we show the posterior of bin_5; as one can notice, our mass modeling is capable of retrieving the true θ space within a 1σ uncertainty. For all velocity bins tested (but not shown here), similarly good results were found. Therefore, we conclude that our mass modeling approach has proven robust to disentangle the various components of RCs.

Fig. B.2. Example of MCMC posteriors of mass-modeled synthetic coadded RCs, in the case of the Burkert and the NFW halo, on the left and right, respectively. Here, we show that the posteriors correspond to the synthetic coadded RC shown in Figure B.1. The best-fit parameters are printed above each posterior with a 1σ uncertainty. The vertical black dashed lines show the best-fit central value along with the 1σ upper and lower limit, and the true values are shown by blue crosses.

Appendix C: Tables

In Table C.1, we provide the physical parameters of observed galaxies (in the case of coadded RCs). In particular, we provide relevant details on coadded RCs (binning), averaged physical quantities (${\stackrel{\sim }{R}}_{\mathrm{D}}$, ${\stackrel{\sim }{V}}_c$, ${\stackrel{\sim }{M}}_{\ast }$, $\stackrel{\sim }{\mathrm{SFR}}$, and ${\stackrel{\sim }{R}}_{\mathrm{H}2}$) per bin, and best-fit modeled parameter space θ and the DM fraction, in the case of the Burkert and NFW halo. In Table C.2, we provide the relevant physical parameters (${\stackrel{\sim }{R}}_{\mathrm{D}}$, ${\stackrel{\sim }{V}}_c$, ${\stackrel{\sim }{M}}_{\ast }$, ${\stackrel{\sim }{R}}_{\mathit{gas}}$, and ${\stackrel{\sim }{f}}_{\mathrm{DM}}$) of galaxy simulations: EAGLE, TNG100, and TNG50.

Appendix D: Extra Figures

In Figures D.1 & D.2, we provide the posterior distribution of full parameter space θ for individual and coadded RCs, respectively.

Fig. D.1. Individual RC mass modeling: A few examples of full posterior distributions (MCMC output) of the estimated parameters for the 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 the case of a Gaussian prior) of parameter space θ. This figure is an extended version of Figure 2.

Fig. D.2. Coadded RC mass modeling: A few examples of full posterior distributions (MCMC output) of the estimated parameters for the 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 the case of the Gaussian prior) of parameter space θ. This figure is an extended version of Figure 3.

All Tables

All Figures

Fig. 1. Coadded and binned rotation curves. The individual and binned RCs are shown by blue thin curves and squares connected with a dark blue line, respectively. The central velocity (${\stackrel{\sim }{V}}_{\mathrm{c}}$) and the number of data points per bin are printed in the legend of each panel. For example, the notation ${\stackrel{\sim }{V}}_{\mathrm{c}}:82$ (15; 64) means that the central circular velocity of the bin is 82 km s−1 and that it contains 15 RCs and a total of 64 data-points. We note that here ${\stackrel{\sim }{V}}_{\mathrm{c}}$ is the circular velocity computed at ${\stackrel{\sim }{R}}_{\mathrm{out}}$ (=2.95 Re).

In the text

Fig. 2. Individual RC mass modeling. Row 1: posterior distributions (MCMC output) of the following estimated parameters: r0 and ρ0 for Burkert (red) and rs and ρs for the NFW (blue) halo. The vertical and horizontal gray dashed lines show the initial guess for parameters r0/s and ρ0/s. The full posterior distribution of parameter space θ of these example cases are shown in Fig. D.1. Rows 2 and 3: best fit to the RCs and their velocity decomposition for Burkert and NFW halos, respectively. In the 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σ 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 the gray shaded area represents the unresolved region of the RC. The best-fit value of parameter space θ is given in the legend of the figures. All the parameters are printed in log, having physical units as follows: MD [M⊙], RD [kpc], Mbulge [M⊙], r0/s [kpc], and ρ0/s [gm cm−3].

In the text

Fig. 3. Coadded RC mass modeling. Row 1: posterior distributions (MCMC output) of the following estimated parameters: ${\stackrel{\sim }{r}}_0$ and ${\stackrel{\sim }{\rho }}_0$ for Burkert (red) and ${\stackrel{\sim }{r}}_s$ and ${\stackrel{\sim }{\rho }}_s$ for the NFW (blue) halo. The vertical and horizontal gray dashed lines show the initial guess for the parameters ${\stackrel{\sim }{r}}_{0/s}$, ${\stackrel{\sim }{\rho }}_{0/s}$, and ${\stackrel{\sim }{M}}_{\mathrm{bulge}}$. The full posterior distribution of parameter space θ of these example cases are shown in Fig. D.2. Rows 2 and 3: best fit to the RCs and their velocity decomposition for the Burkert and NFW halos, respectively. In the 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σ 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 the gray shaded area represents the unresolved region of the RC. The best-fit value of parameter space θ is given in the legend of the figures. All the parameters are printed in log, having physical units as follows: ${\stackrel{\sim }{M}}_{\mathrm{D}}[M_{\odot }]$, ${\stackrel{\sim }{R}}_{\mathrm{D}}[\mathrm{kpc}]$, ${\stackrel{\sim }{M}}_{\mathrm{bulge}}[{\mathrm{M}}_{\odot }]$, ${\stackrel{\sim }{r}}_{0/s}[\mathrm{kpc}]$, and ${\stackrel{\sim }{\rho }}_{0/s}[\mathrm{gm}{\mathrm{cm}}^{-3}]$.

In the text

Fig. 4. Dynamically derived integrated quantities. Left panel: comparison of the photometric and dynamical stellar masses in the Burkert (blue) and NFW (pink) halo profiles for coadded RCs. The black solid line shows the one-to-one relations. We note that the dynamical stellar masses obtained assuming 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 dex. Right panel: dynamically derived bulge masses versus stellar masses, in the case of Burkert (blue) and NFW (pink) halo profiles for coadded RCs. We compared our results of the same redshift range (z = 0.6 − 1.2) and a similar stellar mass sample with a pioneering study by Genzel et al. (2020, red filled circles), and we found that in our sample B/T is never higher than ∼0.3, which is different from what is reported by Genzel et al. (2020).

In the text

Fig. 5. DM fraction within Rout for the Burkert and NFW halos, on the left and right, respectively. The color coding in both panels are the same for the individual RCs (blue filled circles), coadded RCs (blue hexagons), local SFGs (coral red circles), and high-z SFGs (coral red hexagons). The size of the makers in local SFGs includes their errors. The blue shaded area represents the model-independent DM fraction study that we previously conducted in Sharma et al. (2021b).

In the text

	Fig. 6. Concentration vs. virial mass relation for individual and coadded RCs (or galaxies) in blue circles and hexagons, respectively. The blue and red, solid and dashed lines represent the c − Mvir relation of DM-only simulations (Bullock et al. 2001) and the NIHAO suite of hydrodynamical cosmological simulations (Wang et al. 2015; Freundlich et al. 2020), for z = 0 and z ∼ 1, respectively.
In the text

Fig. 7. DM core density vs. core radius relation of local and high-z galaxies. The blue circles and hexagons represent the z ∼ 1 individual and coadded RCs, respectively. The coral red circles and stars represent the local star-forming and low surface brightness galaxies (Persic et al. 1996; Di Paolo et al. 2019), respectively. The brown dashed line is the best fit ${\rho }_0\propto r_0^{\alpha }$ to z ∼ 1 data with a slope α = −1.54, and the orange shaded area represents the 1 − 3σ intrinsic scatter in the relation.

In the text

Fig. 8. DM fraction as a function of the radius in various velocity bins of observed and simulated galaxies. The central value of fDM(< 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σ error in each measurement is represented by a shaded area; whereas in the observations, it is a dashed line following the same color code. The gray shaded area in each panel represents the unresolved region in the observations. The horizontal red dotted-dashed line in each panel shows the 51% DM fraction. The central value of the circular velocity of each bin is printed at the top of each panel, where ${\stackrel{\sim }{V}}_c^{\mathrm{SIM}}$ and ${\stackrel{\sim }{V}}_c^{\mathrm{OBS}}$ specify the simulations and observations, respectively.

In the text

Fig. 9. Structural properties of the DM halo at z ∼ 1. Top: DM core radius and core density as a function of the circular velocity shown in the upper and lower panels, respectively. The individual and coadded RC measurements are represented by small blue circles and hexagons, respectively. The comparison sample of local LSBs (Di Paolo et al. 2019) and SFGs (Persic et al. 1996; Lapi et al. 2018) 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 r0 − Vc plane and below in the ρ0 − Vc plane with respect to the z ∼ 1 SFGs. This is due to the fact that DM cores of z ∼ 1 SFGs are smaller and denser than their local counterparts, indicating an evolution in the structure of DM halos over cosmic time. Bottom: ratio of DM core radius and density at z ≈ 0 to z ∼ 1 as a function of circular velocity. This plot shows the expansion factor of DM cores at z ∼ 1 with respect to local SFGs. The color codes for data points are similar to those in the top 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 a solid line shows the 1σ intrinsic scatter in the relations.

In the text

Fig. A.1. Comparison of the physical parameters of the 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, the so-called main sequence of SFGs. Middle panel: Effective radii of galaxies as a function of stellar mass, the so-called mass-size relation. Right panel: Distribution of the maximum circular velocity of the simulated sample compared to Vout = Vc of the observations.

In the text

	Fig. A.2. Individual 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 the velocity bin is shown in the title of each plot.
In the text

	Fig. A.3. Coadded and binned circular velocity curves of galaxy simulations from first to last EAGLE, TNG100, and TNG50. The color for each panel (and plot) is given in the legend as well as the name of the simulation, and the velocity bin is shown in the title of each plot.
In the text

	Fig. A.4. Binned RCs of all velocity bins of galaxy simulations from first to last EAGLE, TNG100, and TNG50. The color maps give an idea of their increasing circular velocities.
In the text

Fig. B.1. Example of synthetic coadded RCs. Top panel: Observed coadded and binned RC, which is shown as a reference. Bottom panel: Coadded and binned synthetic RCs in the case of the Burkert and the NFW dark matter halo, which are on the left and right, respectively. The color code in all of the panels are the same and given as follows: the brown square represents the individual RCs and the blue stars represent the coadded RC.

In the text

Fig. B.2. Example of MCMC posteriors of mass-modeled synthetic coadded RCs, in the case of the Burkert and the NFW halo, on the left and right, respectively. Here, we show that the posteriors correspond to the synthetic coadded RC shown in Figure B.1. The best-fit parameters are printed above each posterior with a 1σ uncertainty. The vertical black dashed lines show the best-fit central value along with the 1σ upper and lower limit, and the true values are shown by blue crosses.

In the text

	Fig. D.1. Individual RC mass modeling: A few examples of full posterior distributions (MCMC output) of the estimated parameters for the 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 the case of a Gaussian prior) of parameter space θ. This figure is an extended version of Figure 2.
In the text

	Fig. D.2. Coadded RC mass modeling: A few examples of full posterior distributions (MCMC output) of the estimated parameters for the 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 the case of the Gaussian prior) of parameter space θ. This figure is an extended version of Figure 3.
In the text