Rotation Curves in z ∼ 1–2 Star-forming Disks: Evidence for Cored Dark Matter Distributions

Stellar and gas (H i, CO, Hα) rotation curves (RCs) have been a valuable tool for investigating the structure and mass distribution of galaxies, ever since the first stellar kinematic studies of the Milky Way (MW) in the early part of the last century (Kapteyn 1922; Oort 1927, see Bland-Hawthorn & Gerhard 2016). This continued in the 1960–2000s in disk galaxies in the nearby universe (Freeman 1970; Rubin & Ford 1970; Rogstad & Shostak 1972; Roberts & Rots 1973; Roberts & Whitehurst 1975; Rubin et al. 1978; van der Kruit & Allen 1978; Persic & Salucci 1988; Begeman et al. 1991; Sofue & Rubin 2001; Catinella et al. 2006; de Blok et al. 2008; Lelli et al. 2016). RCs with flat or positive slopes (circular velocity v_c increasing with radius R) on scales of tens of kiloparsecs have since been one of the fundamental pillars of the dark matter paradigm, assuming that Newtonian physics applies on large scales (Einasto et al. 1974; Ostriker et al. 1974; Kent 1986; Courteau et al. 2014).

The overall RC (baryons plus dark matter) is determined by the mass fraction of the baryons in the disk and bulge, to dark matter in the halo, m_b = (M_bulge + M_disk,* + M_disk,gas)/M_DM, by the concentration of the halo, and by the specific angular momentum of the halo and the baryons. The angular momentum parameter of the baryons at the disk scale (Appendix A.2), λ_baryon, is often assumed to be similar to that of the dark matter at the virial scale (λ_DM). This simple assumption is supported by observations (Romanowsky & Fall 2012; Fall & Romanowsky 2013; Burkert et al. 2016), but theoretically not necessarily expected (Übler et al. 2014; Danovich et al. 2015; Teklu et al. 2015; Jiang et al. 2019). If the baryons are in a thin disk with an exponential distribution (Sérsic index n_S = 1), the half-mass (or effective) radius of the exponential disk is ${R}_{e}\sim ({\text{}}1.68/\sqrt{2})\times {\lambda }_{\mathrm{baryon}}\times {R}_{\mathrm{virial}}$ (Mo et al. 1998).

For the MW, the observed RC rises from the center to a peak at 6 kpc, then exhibits a shallow decline to 25 kpc, and then is flat from 25 to 60 kpc (left panel of Figure A1; Battacharjee et al. 2014; Bland-Hawthorn & Gerhard 2016; Eilers et al. 2019; Reid et al. 2019; Sofue 2020). The MW ratio of dark to total mass at R_e ∼ 4.5 kpc is about f_DM(R_e) = 0.38 ± 0.1 (Bovy & Rix 2013; Bland-Hawthorn & Gerhard 2016) such that the MW disk is baryon dominated at R_e and becomes dark matter dominated at the solar circle. The uncertainties of the various estimates are large enough not to rule out that the MW is a "maximal" disk (f_DM(R_e) ≡ 0.28; van Albada et al. 1985). The MW RC is typical of other massive, bulged disks at z ∼ 0 (Sofue & Rubin 2001; van der Kruit & Freeman 2011). Very massive (v_c ∼ 250–310 km s⁻¹) disks with large bulge to total ratios have f_DM(R_e) = 0.25 ± 0.15 but are rare (Barnabe et al. 2012; Dutton et al. 2013). Otherwise, typical disk galaxies in the local universe are more dark matter dominated, with f_DM(R_e) ∼ 0.5–0.9. The most dark-matter-dominated disks tend to have low baryonic mass and circular velocity (Martinsson et al. 2013a, 2013b; Courteau & Dutton 2015).

The cosmic star formation density peaked 5–11 Gyr ago (z ∼ 1–2.5; Madau & Dickinson 2014), during which galaxy halos containing MW mass galaxies first formed in large numbers (e.g., Mo & White 2002). Over the past two decades, high-throughput, adaptive-optics-assisted, near-integral field spectrometers (IFS), such as SINFONI on the ESO-VLT (Eisenhauer et al. 2003; Bonnet et al. 2004) or OSIRIS on the Keck telescope (Larkin et al. 2003), and seeing-limited, multiplexed IFSs, such as KMOS at the VLT (Sharples et al. 2012), have become available on 8–10 m telescopes. With these IFSs, it is now possible to carry out deep, velocity-resolved (FWHM ∼ 80–120 km s⁻¹) spectroscopic imaging of Hα in z ∼ 0.6–2.6 disk galaxies on the star-forming "main sequence (MS)" (Whitaker et al. 2012, 2014; Speagle et al. 2014).

Over the past decade, we have carried out two main IFS surveys of high-z galaxy kinematics, SINS and zC-SINF with SINFONI (Förster Schreiber et al. 2006, 2009, 2018; Genzel et al. 2006) and KMOS^3D with KMOS (Wisnioski et al. 2015, 2019). In total, we have assembled ∼850 IFS data sets of MS star-forming galaxies (SFGs) across the cosmic star/galaxy formation peak, covering the mass range from log(M_∗/M_⊙) = 9.5–11.5. Between 60% and 80% of the MS galaxies at these masses are rotation-dominated, turbulent, and thick disks with v_rot(R_e)/σ₀ ∼ 3–10 (Wisnioski et al. 2015, 2019; Simons et al. 2017). Here, v_rot(R_e) is the inclination- and beam-smearing-corrected intrinsic rotation velocity of the disk at R_e, and σ₀ is the average velocity dispersion of the (outer parts of the) disk, after removal of beam-smeared rotation, instrumental line broadening, and other orbital motions. In our SINS and KMOS^3D surveys, we have emphasized deep integrations, for high-quality data on individual galaxies, rather than optimizing the number of observed galaxies—an approach that has enabled excellent-quality high-z RCs.

In addition RCs of high-z galaxies can now also be obtained from millimeter and submillimeter emission lines, such as CO and [C ii] (e.g., Genzel et al. 2013; Übler et al. 2018). With the availability of such high-quality high-z RCs, it is possible to go beyond the mere demonstration of rotationally dominated motions (Förster Schreiber & Wuyts 2020) to exploiting RCs for studying mass distributions in more detail.

Wuyts et al. (2016) modeled the inner disk kinematics with the first and second velocity moments of Hα along the kinematic major axis in 240 z = 0.6–2.6 SFGs from the KMOS^3D survey (Wisnioski et al. 2015, 2019). They compared the dynamical masses within R_e with the sum of stellar and gas masses inferred from multiband optical/IR photometry, Hubble Space Telescope (HST) imaging, and molecular gas scaling relations (e.g., Tacconi et al. 2018). Wuyts et al. (2016) found average baryon fractions of f_baryon(<R_e) ∼ 60%, increasing with redshift to ∼90% at z ∼ 2.2. Übler et al. (2017) inferred a similar global increase in baryon fractions from z < 1 to z > 2 for 135 "best disks" in a "Tully–Fisher offset analysis (TFA)" of KMOS^3D SFGs. Turner et al. (2017) carried out a TFA study in 1200 z = 0–3 SFGs in the KMOS-KDS and literature samples, and Tiley et al. (2016, 2019a) in 210 z = 1 SFGs in KMOS-KROSS. All three TFA studies find similar offsets as a function of z. Van Dokkum et al. (2015) observed integrated Hα line widths in 25 compact SFGs with MOSFIRE at the Keck telescope and concluded that these widths are consistent with a Keplerian falloff of circular velocities, indicative of masses dominated by dense stellar components.

Genzel et al. (2017) presented deep SINFONI and KMOS^3D Hα imaging spectroscopy for six very massive, bulgy z = 0.9–2.4 SFGs (log(M_baryon/M_⊙) = 11.1–11.3), tracing individual RCs to 1.5–3 R_e for the first time at these redshifts. In four or five of the six SFGs, the RCs decline beyond a maximum near R_e. Genzel et al. (2017) interpreted these radial drops as being due to a combination of asymmetric drift at high velocity dispersions and large baryon fractions within R_e (f_DM(R_e) = 0 $\to$ 30%).

To test whether the results of Genzel et al. (2017) are representative of the overall population of high-z massive SFGs, Lang et al. (2017) stacked 101 massive, large, rotationally supported SINFONI and KMOS^3D SFGs with median z ≥ 1.52 and log(M_∗/M_⊙) = 10.6, which had robustly detected turnover radii. They normalized individual RCs by turnover radius and corresponding velocity before stacking. The resulting average RC shows a radially decreasing RC, in agreement with Genzel et al. (2017). Lang et al. (2017) also showed that a second method of "stellar light normalization" using the effective radius R_e and Sérsic index n from HST NIR imaging to estimate the turnover radius gave comparable results. Using four different approaches to estimate dark matter fractions near R_e, the results of Genzel et al. (2017), Lang et al. (2017), Übler et al. (2017), and Wuyts et al. (2016) all indicate that massive rotating disks have low dark matter fractions within several R_e at z > 1–2.5. Dark matter fractions increase toward lower redshift, in agreement with the local universe results described above.

Other Studies. Price et al. (2016, 2020) analyzed Hα slit kinematics from the Keck MOSDEF survey for 681 z ∼ 1.4–3.8 SFGs, for another independent look at the inner disk kinematics. Their results are in excellent agreement with Wuyts et al. (2016). Tiley et al. (2019b) compiled 1500 z = 0.6–2.2 SFGs from the KROSS (551 galaxies with median z = 0.85 and log(M_∗/M_⊙) = 10.0), KGES (228 galaxies with median z = 1.5 and log(M_∗/M_⊙) = 10.3), KMOS^3D (145 galaxies with median z = 2.3 and log(M_∗/M_⊙) = 10.3) surveys, as well as MUSE observations (96 galaxies with median z = 0.67 and log(M_*/M_⊙) = 9.8). Tiley et al. confirmed dropping RCs similar to Lang et al. (2017) for all their substacks when adopting a similar "self-normalization" methodology based on the RC turnover but found instead flat or even rising average RCs when normalizing by the observed velocity at six times the disk scale length, 6R_d, assuming exponential light distributions (and accounting for spatial beam smearing). This illustrates the impact of different approaches, and possibly sample differences, in stacking and highlights the need to study individual RCs.

In this paper, we present 41 high-quality, individual z = 0.65–2.45 RCs (35 new RCs and the 6 RCs from Genzel et al. 2017) that we have collected and analyzed since Genzel et al. (2017), which we will refer to as the "RC41" sample. In following up on Genzel et al. (2017), our main goal was to increase statistics and to expand the sample to lower masses and lower redshifts. With RC41, we also better sample the overall MS population across the peak of cosmic galaxy formation than in the earlier work. Of the 41 RCs, 17 come from SINFONI at the ESO-VLT in both non-AO and AO (14) modes. For seven, we combined SINFONI observations with data obtained with the KMOS multi-IFS instrument at the VLT. For eight, we used KMOS data only. We also compare and combine Hα- and molecular gas (CO-based) RCs in two galaxies. Genzel et al. (2013) and Übler et al. (2018) have shown excellent agreement of CO (from IRAM-NOEMA) and Hα (from LBT-LUCI) RCs for two SFGs with comparable spatial resolutions. In RC41, we now add seven new NOEMA-CO RCs. RC41 on-source integration times varied from 4 to 56 hr, with a median of 16 hr, and a total of 800 on-source hours, equivalent to at least 1000 total telescope hours on the ESO-VLT, on the LBT, and on NOEMA. Figure 1 shows the locations of the RC41 sources in the planes of log M_∗ versus MS offset (δMS = log(SFR/SFR(MS, z))), effective radius R_e of optical (5000 Å) stellar continuum, baryonic (sum of stellar and gas) mass surface density within R_e, log Σ_baryon, and bulge mass M_bulge. Throughout this paper, filled blue and red circles denote RC41 galaxies in the redshift slices 0.65–1.2, and 1.2–2.45. Table 1 summarizes the relevant basic "input" properties of the RC41 sample.

Zoom In Zoom Out Reset image size

Table 1.  Overview of the RC41 Sample and Its "Input" and Environmental Properties

	Galaxy	z	Morphology	Inclination	vc(Re(A)) km s−1	Re (A) kpc	kpc arcsec−1	logM* HST Input (M⊙)	B/T Optical Light	B/T Dynamical	Satellite/Companion/Merger	Distance (arcsec)	Light/Mass Ratio Central/Satellite	Δv (km s−1) Satellite–Central	δv/vA in Interaction	RJ/RAB
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)
1	EGS3_10098/ EGS4_30084	0.66	outer spiral/ring+bulge	31	326	3	6.9	11.1	⋯	0.60	⋯	⋯	⋯	⋯	⋯	⋯
2	U3_21388	0.67	N3079 disk+bulge	82	230	7.0	7.0	10.8	0.39	0.05	⋯	⋯	⋯	⋯	⋯	⋯
3	EGS4_21351	0.73	bulge and ring	47	187	3.3	7.2	10.9	0.21	0.46	⋯	⋯	⋯	⋯	⋯	⋯
4	EGS4_11261	0.75	bulge & clumpy disk	60	327	4.0	7.3	11.3	⋯	0.50	⋯	⋯	⋯	⋯	⋯	⋯
5	GS4_13143	0.76	spiral +clumpy nuclear disk	74	156	5.6	7.4	9.8	⋯	0.70	interaction	2'' ENE	H: 2.2:1	165	0.11	0.15
6	U3_05138	0.81	bulge and disk	55	145	7.5	7.5	10.2	0.03	0.50	companion	13 NNW	H: 50:1	⋯	⋯	⋯
7	GS4_03228	0.82	bulge and disk	78	146	7.1	7.6	9.5	⋯	0.80	⋯	⋯	⋯	⋯	⋯	⋯
8	GS4_32976	0.83	bulge and disk	68	249	6.8	7.6	10.4	0.09	0.90	⋯	⋯	⋯	⋯	⋯	⋯
9	COS4_01351	0.85	bulge+ring	68	290	8.0	7.7	10.7	0.02	0.24	⋯	⋯	⋯	⋯	⋯	⋯
10	COS3_22796	0.91	small bulge & big ring	58	145	9.0	7.8	10.3	0.1	0.15	⋯	⋯	⋯	⋯	⋯	⋯
11	U3_15226	0.92	bar-wide spiral	50	194	5.5	7.9	11.1	0.24	0.55	⋯	⋯	⋯	⋯	⋯	⋯
12	GS4_05881	0.99	clumpy disk +extremely red bulge?	60	200	5.6	8.0	9.8	⋯	0.85	⋯	⋯	⋯	⋯	⋯	⋯
13	COS3_16954	1.03	bulge+spiral arms	49.5	282	8.1	8.1	10.7	0.08	0.50	companion?	3'' WNW	H: 26:1	⋯	⋯	⋯
14	COS3_04796	1.03	spiral+bulge	50	269	9.7	8.1	10.8	0.04	0.18	satellite	11 ENE	H 13:1	200	0.08	0.25
15	EGS_13035123	1.12	spiral+bulge	24	220	10.2	8.2	11.2	0.06	0.20	⋯	⋯	⋯	⋯	⋯	⋯
16	EGS_13004291	1.20	clumpy disk,+ large bulge, strong tidal arm	27	354	3.0	8.3	11.0	⋯	0.61	late stage merger	⋯	⋯	⋯	⋯	⋯
17	EGS_13003805	1.23	spiral+big extincted bulge	37	393	5.6	8.3	11.2	⋯	0.29	⋯	⋯	⋯	⋯	⋯	⋯
18	G4-38153	1.36	edge on ring	75	356	5.9	8.4	10.4	⋯	0.16	⋯	⋯	⋯	⋯	⋯	⋯
19	G4_24985	1.40	extincted bulge and disk	40	300	4.6	8.4	10.9	⋯	0.40	⋯	⋯	⋯	⋯	⋯	⋯
20	zC_403741	1.45	bulge and ring	28	204	2.6	8.5	10.7	⋯	0.68	⋯	⋯	⋯	⋯	⋯	⋯
21	D3a_6397	1.50	bulge and ring	30	308	6.4	8.5	11.1	⋯	0.57	⋯	⋯	⋯	⋯	⋯	⋯
22	EGS_13011166	1.53	small bulge & clumpy spiral arms	60	348	6.3	8.6	11.1	⋯	0.55	⋯	⋯	⋯	⋯	⋯	⋯
23	GS4_43501 (GK 2438)	1.61	bulge and ring	62	257	4.9	8.5	10.6	0.06	0.40	⋯	⋯	⋯	⋯	⋯	⋯
24	GS4_14152	1.62	bulge and disk	55	397	6.8	8.5	11.3	⋯	0.23	⋯	⋯	⋯	⋯	⋯	⋯
25	K20_ID9-GS4_27404	2.04	small bulge and ring	48	250	7.1	8.3	10.6	0.14	0.30	⋯	⋯	⋯	⋯	⋯	⋯
26	zC_405501	2.15	bulge & edge on clumpy disk	75	139	5.0	8.3	9.9	⋯	0.07	interacting group with 3 members	N+17 SSE −175	H: 2.7:1 5.4:1	⋯	⋯	⋯
27	SSA22_MD41	2.17	ring	72	189	7.1	8.3	9.9	⋯	0.05	⋯	⋯	⋯	⋯	⋯	⋯
28	BX389	2.18	edge on ring	76	331	7.4	8.3	10.6	⋯	0.30	satellite	07 SSW	m: 13:1	180	0.12	0.27
29	zC_407302	2.18	disk	60	280	4.0	8.3	10.4	⋯	0.50	satellite	055 N	m: 10:1	90	0.19	0.22
30	GS3_24273	2.19	extincted bulge & clumpy disk	60	222	7.0	8.3	11.0	0.15	0.80	⋯	⋯	⋯	⋯	⋯	⋯
31	zC_406690	2.20	empty clumpy ring	25	301	4.5	8.3	10.6	⋯	0.90	interacting fluffy satellite	16 W	m: 8:1	100?	⋯	⋯
32	BX610	2.21	extincted bulge and disk	39	327	4.9	8.2	11	⋯	0.42	⋯	⋯	⋯	⋯	⋯	⋯
33	K20_ID7-GS4_29868	2.23	empty clumpy ring	64	308	8.2	8.2	10.3	⋯	0.03	two companions?	32 S 31 W	H: 4:1	⋯	⋯	⋯
34	K20_ID6-GS3_22466-GS4_33689	2.24	disk	31	162	5.0	8.2	10.4	⋯	0.30	minor merger	13 SW	H: 6:1	−170	0.05	0.15
35	zC_400569	2.24	bulge and disk	45	289	4.0	8.2	11.1	⋯	0.70	group with 3 members	12 S 25 S	m: 11:1	−280 and −400	0.01	0.05
36	BX482	2.26	small extincted bulge and ring	60	293	5.8	8.2	10.3	⋯	0.02	interacting group with 3 members	19 SW 33 SE	H: 10:1 20:1	690 and 800	0.01	0.09
37	COS4_02672	2.31	bulge and ring	62	190	7.4	8.2	10.6	⋯	0.10	companion	26 ENE	H 5:1	⋯	⋯	⋯
38	D3a_15504	2.38	bulge and disk	40	266	6.1	8.1	11.0	⋯	0.30	interacting satellite	15 NW	m: 30:1	−50	0.06	0.08
39	D3a_6004	2.39	big bulge and ring	20	416	5.3	8.1	11.5	⋯	0.44	companion?	17 SE	H: 50:1	⋯	⋯	⋯
40	GS4_37124	2.43	bulge and red ring/disk	67	260	3.2	8.1	10.6	0.20	0.70	⋯	⋯	⋯	⋯	⋯	⋯
41	GS4_42930 (GK 2363)	2.45	bulge and compact disk	59	171	2.8	8.1	10.3	0.25	0.50	⋯	⋯	⋯	⋯	⋯	⋯

Note. We are adopting a Ω_m = 0.3, H₀ = 70 km s⁻¹ Mpc⁻¹ ΛCDM universe and a Chabrier (2003) initial stellar mass function. Effective radii, inclinations, and input stellar masses are derived from analysis of the optical HST ancillary information, using the techniques of Wuyts et al. (2016) (see the text). Column 15: δv/v_c(R_e(A)) = (v_c(R_e(A))/Δv) × (R_e(A)/R_AB) × (M_B/M_A); column 16: R_J/R_AB(deprojected) = (R_e(A)/1.5 × R_AB) × (M_B/2M_A)^0.333.

Download table as:  ASCIITypeset images: 1 2

Our VLT targets are primarily selected from and benchmarked to the 3D-HST and CANDELS surveys (Grogin et al. 2011; Koekemoer et al. 2011; Skelton et al. 2014; Momcheva et al. 2016). Additional targets come from the Steidel et al. (2004) "BX-BM" Lyman-break surveys, the GMASS Spitzer survey (Cimatti et al. 2008), the zCOSMOS-Deep survey of the COSMOS field (Lilly et al. 2007; Scoville et al. 2007), and BzK galaxies from the Deep3a field (Kong et al. 2006). The observations were taken as part of the SINS/zC-SINF Hα survey of 92 galaxies with SINFONI (Förster Schreiber et al. 2009, 2018; Mancini et al. 2011) and the KMOS^3D survey of 740 galaxies with KMOS (Wisnioski et al. 2015, 2019), each with an Hα detection fraction of about 80% (and very small overlap of 18 sources observed in both surveys). At NOEMA in the northern hemisphere, the targets were observed in the PHIBSS1 & 2 CO surveys (Tacconi et al. 2010, 2013, 2018; Übler et al. 2018), selected and benchmarked from the All-Wavelength Extended Groth Strip International Survey (abbreviated here as EGS: Davis et al. 2007; Noeske et al. 2007), with spectroscopic redshifts from DEEP2 (Newman et al. 2012). Additional PHIBSS targets were taken from the zCOSMOS-Deep survey (Lilly et al. 2007) and the rest from 3D-HST in the GOODS-N, EGS, and COSMOS fields. Most of the 180 PHIBSS targets were observed in compact array configurations as the focus was on source-integrated properties, with an 82% detection fraction; 12 of the CO-detected galaxies were followed up in extended configurations to more fully resolve their line emission.

Sample Selection. From the full Hα and CO parent sample of 800 unique detected galaxies, with one exception, we selected rotationally supported (v_rot/σ₀ > 2.3) SFGs with stellar masses 9.5 ≤ log(M_*/M_⊙) ≤ 11.5, redshifts 0.65 ≤ z ≤ 2.5, and near the MS, −0.6 ≤ δMS ≤ 1. Given the atmospheric transparency in the near-infrared and millimeter ranges, this splits the data into three redshift slices, 0.65 ≤ z ≤ 1.2, 1.2 ≤ z ≤ 1.6, and 2.02 ≤ z ≤ 2.5. These selections ensure a fairly homogeneous coverage of the star-forming disk population around the MS in the stellar mass–star formation rate plane in each redshift slice. The kinematic criterion removes about one-fourth of the detected galaxies; of the remaining objects, 86% (514) satisfy the stellar mass, redshift, and MS offset cuts.

We retained only spatially well-resolved data sets, necessary to derive RCs. This restriction culls compact SFGs (R_e < 2 kpc), given the FWHM angular resolution of our data sets, which range from 025 for AO data sets to 08 for the lowest-resolution KMOS and NOEMA data, leaving 415 objects. In addition, we required significant line detection at R ≥ R_e and further eliminated sources for which residual OH sky-line emission affected the Hα line profile over part of the galaxy, leading to the final RC41 sample that represents the best cases in terms of signal-to-noise ratio, radial coverage, and quality among our data sets. The statistics are best at the massive tail of the star-forming MS population (bottom left of Figure 1), and the median optical effective radius in RC41 is 5.7 kpc. The one exception noted above is a more actively star-forming disk at δMS ∼ 1.25, which otherwise meets all other criteria applied including the particularly stringent ones of well-detected and well-resolved line emission out to large radii. Taken together, the two size and line emission extent criteria result in the largest cut and bias the RC41 sample to larger disks (bottom-right panel of Figure 1), which is an inevitable but acceptable consequence of wanting to push high-z RCs to ≥0.1 R_virial. This means that the baryonic surface densities of the disks in RC41 are at the low-density tail of the overall Σ_baryon distribution (upper-left panel of Figure 1). With the main conclusions of this work in mind, it is fair to say that this selection is unlikely to bias dark matter fractions downward. As discussed in van Dokkum et al. (2015) and Wuyts et al. (2016), more compact SFGs, missing from RC41, have higher baryonic surface densities and are likely to be more baryon dominated than RC41 galaxies.

Dynamical Modeling. As in our earlier work (Genzel et al. 2006, 2017; Burkert et al. 2016; Wuyts et al. 2016; Lang et al. 2017; Übler et al. 2017, 2018), we use forward modeling from a parameterized, input mass distribution to establish the best-fit models for a given Hα or CO data set. This mass model is the sum of an unresolved passive bulge (not emitting in Hα or CO), a rotating flat disk of Sérsic index n_S, effective (half-light) radius R_e and (constant) isotropic velocity dispersion σ₀, and a surrounding halo of dark matter. As discussed in more detail in the above-cited references and in Appendix A, we compute from this mass three-dimensional data cubes (I(x, y, v_z)) of the disk gas, convolved with a three-dimensional kernel describing the instrumental point-spread function PSF (δx, δy, δv_z) of our measurements. We then compare directly this "beam-smeared" model to the observed data and vary model parameters to obtain best fits. The most common approach (taken in this paper) is to extract velocity centroids (from Gaussian fits) and velocity dispersion cuts from a suitable software slit along the dynamical major axis of the galaxy, for both the model and measurement cubes. We either use a constant software slit width (typically ∼1–1.5 FWHM of the data set), or for low-inclination targets where the isovelocity contours fan out, a fanned slit width (5°–10° opening angle; e.g., van der Kruit & Allen 1978). This 1D method is applicable to all RC41 galaxies, as most of the RC information is contained along this major axis (see Genzel et al. 2006, 2017). For the highest resolution data sets, or very large galaxies, we obtain additional information by investigating the model and measurement 2D velocity and velocity dispersion maps, or even comparing the individual spaxels in the data cubes for a full 3D analysis. We discuss these 2D and 3D analyses for suitable RC41 galaxies in Paper II of this series (S. Price et al. 2020, in preparation) and refer the reader for the full detailed discussion to that paper.

Figure 2 shows three examples of the full data sets for each galaxy, along with the best-fit model (red line in the right three panels), obtained from the analysis described in Appendix A. Figure 3 shows the inclination-corrected model RCs of the baryons, dark matter, and total mass for the three galaxies in Figure 2. For each, we have two- or three-band HST imagery, plus high-quality Hα or CO brightness distributions, and first- and second-moment kinematic maps. The bottom (Figure 2) case (left panel in Figure 3) is a strongly baryon-dominated z ∼ 2 galaxy (also contained in Genzel et al. 2006, 2017), where the overall intrinsic circular velocity curve drops with radius and where turbulent pressure corrections (asymmetric drift, A1) are very important and lead to the precipitous drop in rotation velocity. In contrast, the galaxy shown in the middle panels is a dark-matter-dominated z ∼ 2 system, with a rising RC. The top panels in Figure 2 (right panel in Figure 3) feature a larger, z ∼ 1 galaxy with a flat RC.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The 41 best major axis position–velocity diagrams make up RC41. In this paper, we focus on the analysis of the 1D cuts of the Hα/CO data on the kinematic and structural major axis. Figure 4 summarizes all of the observed RCs of the RC41 sample (black filled circles with rms uncertainties), along with the projected (and beam-convolved) best-fit, model RCs (red curves). Tables D1 and D2 in Appendix D summarize the inferred properties of the RC41 galaxies obtained from the analysis. Paper II will examine the 2D distributions and analyze non-axisymmetric distributions and/or radial motions.

Zoom In Zoom Out Reset image size

Shapes and Symmetry of the Observed Rotation Curves. Figure 4 shows that the RCs of RC41 have a variety of shapes. We make a qualitative assessment of their shapes: about half (17 of 41) are falling or exhibit a clearly defined maximum and begin to drop beyond that. About half of these again are at z > 2.1. Eight of the 17 exhibit substantial drops in rotation velocity (factors of 1.2–2 from the peak to the outermost points at ∼2–3 R_e). One-third (15 of 41) of the RC41 SFGs have flat RCs. RCs rising to the outermost measurement point (typically 2 R_e and in three cases 2.5 R_e) are relatively rare (9 of 41, or 20%) but are seen across the full redshift range. For details of the modeling of the individual RCs, including the velocity dispersion, Hα flux cuts, and HST imagery, we refer the reader to Appendix A.

We determined how many of the 41 RCs show significant deviations from reflection symmetry around the dynamical center,¹⁴ and whether the RCs are flat, rising, or reach a peak and turn down, or even fall substantially (by a factor of >1.1) from the maximum to the outermost point. Deviations from point symmetry could be indications of significant perturbations by nearby satellite galaxies, or of warps in the outer galaxy disk. We find that for 33 of the 41 RCs (82.5%), a parabola fit determined from the RC on one side of the galaxy also fits the other side with χ² ≤ 1.2 (Figure 5). Only three (five) SFGs have χ² > 2 (≥1.5). The reflection symmetry of our RC41 sample is quite striking. Simulated high-z galaxies in the currently highest resolution cosmological hydro-simulations, Illustris-TNG50, are far less symmetric (Pillepich et al. 2019; Übler et al. 2020).

Zoom In Zoom Out Reset image size

Environment and Perturbations. Previous observations and simulations indicate that the rates of galaxy interactions and mergers increase rapidly with redshift (Le Fèvre et al. 2000; Conselice et al. 2003; Wetzel et al. 2009; López-Sanjuan et al. 2010; Rodriguez-Gomez et al. 2015; Mantha et al. 2018; Cibinel et al. 2019; Pillepich et al. 2019). For instance, the most recent TNG50 simulations find that at z > 2, even moderately massive galaxies are asymmetric with disturbed disk rotation (Pillepich et al. 2019; Übler et al. 2020). What can we say about the environment and the level of perturbations for the RC41 sample?

Table 1 summarizes our findings. In RC41, there are only two major mergers/interacting systems (mass ratio <3:1), one early stage (GS4_13143, z = 0.76, separation 2'', or 15 kpc), and one late stage, with a prominent tidal tail surrounding a compact, but apparently well settled, compact bulge-disk, EGS_13004291 (z = 1.197, Tacconi et al. 2013). There are eight systems (seven at z > 2.1) where the Hα velocities show evidence for a clear dynamical interaction with a smaller mass galaxy—in those cases, the galaxy mass ratios range from 6:1 to 30:1. Most of these are probably satellites or minor mergers; there are three galaxy groups with two or three members. There are five additional cases of a central galaxy with a nearby projected companion (mass ratio 5:1–50:1) but no spectroscopic evidence for physical association. In summary, 10–15 out of the 41 systems show evidence for companions or satellites. For comparison, an estimated 5%–10% of all KMOS^3D galaxies (z = 0.6–2.6; N. M. Förster Schreiber 2020, private communication) have interacting satellites/minor mergers with mass ratios <30:1. This is consistent with RC41 because the RC41 data are deeper, thus making physical interactions easier to detect. RC41 also has relatively more z > 2 galaxies than the full KMOS^3D sample.

Such interactions and their resulting tidal forces can perturb the motions in the central galaxy and/or induce warps. The last two columns in Table 1 provide approximate estimates of the impact of the satellites/minor mergers on the gas motions at ≥R_e in the central galaxies. In column 15, we compute the relative change in velocity of a particle in central galaxy A, normalized to the circular velocity at R_e, caused by the gravitational pull from satellite B during the "encounter time" δt ∼ 1–1.5 × R_AB/Δv. Here, R_AB is the projected separation between A and B, and Δv is the observed encounter velocity from column 14. The estimate in column 15 is a lower limit for the effect of B on A during a one-time encounter on an unbound or barely bound interaction. In column 16, we compute the ratio of R_e(A) to the "tidal," or Jacobi (or Hills), radius of a particle in galaxy A, which is strongly perturbed by B, given the separation R_AB. If this ratio exceeds unity, then a particle in the central galaxy at R_e will be strongly perturbed in its long-term orbit due to tidal forces from orbiting B (Binney & Tremaine 2008, chapter 8, Equation (8.14)). This estimate gives an upper limit to the effect of tidal forces of the minor merger or satellite, as it assumes a stable, near-circular orbit. Inspection of columns 15 and 16 show that given the mass ratios and separations of the interaction partners, the tidal perturbations in RC41 are modest (∼0.1) and could on average explain the ∼20 km s⁻¹ perturbations in the RCs. This can be significant in some galaxies but is comparable to the measurement uncertainties in the outer disks of most of our galaxies. In further support of this conclusion, of the 10–11 galaxies with spectroscopic evidence for interactions, only 2 have somewhat asymmetric RCs, with χ² > 2 (Figure 5). We cautiously conclude that for most of the RC41 sample there is no convincing evidence that perturbations and environmental interactions affect the kinematics.

Warps. Significant m = 1, 2, and 3 mode, out-of-plane warps (δz ∼ sin(m(φ–φ₀)) are frequently observed in the extended outer H i layers of z ∼ 0 disks, including the MW (van der Kruit & Freeman 2011). The MW H i warp starts at R ∼ 15 kpc, or 4.4 R_e(MW) (Bland-Hawthorn & Gerhard 2016), is a superposition of m = 0, 1, and 2 modes and reaches an average (maximum) scale height of 2 (3) kpc at R = 30 kpc (Levine et al. 2006). Warps in H i may be fire-hose or buckling instabilities near and outside the (possible) truncation of the stellar disk (Binney & Tremaine 2008, chapter 6, van der Kruit & Freeman 2011). Disk stability analysis shows that such instabilities occur below a critical wavelength of λ < λ_crit = σ²/GΣ_disk, where σ is the one-dimensional, in-plane velocity dispersion of the disk and Σ_disk its mass surface density (Toomre 1964; Binney & Tremaine 2008, chapter 6, van der Kruit & Freeman 2011). In hydrostatic equilibrium, the z scale height is given by h_z = σ_z²/GΣ_disk. This means that the buckling/bending instability is only effective in cold disks with σ_z/σ < a_crit ∼ 0.3–0.6 (Toomre 1964; Araki 1985; Merritt & Sellwood 1994). Because z ∼ 1–2.5 SFGs have σ_z/σ ∼ 1 (Wisnioski et al. 2015), the buckling instability should be suppressed or less effective. The reflection symmetry of most RC41 RCs excludes strong even modes (m = 0, 2). Increasing or decreasing rotation velocities on both sides of the rotation curve, with equal probability, can occur if the dominant bending mode is odd (m = 1 or m = 3). Finally, the phase of the warp need not be aligned with the major axis and might change with radius. Such precessions should be detectable in the residual velocity maps. Those galaxies that show such residuals, suggesting warps or radial motions, will be discussed in Paper II.

Summarizing, we find that most of the RCs of the RC41 sample are symmetric in terms of their reflection around the dynamical center. About half of the RCs exhibit a well-defined inner peak and a drop of the rotation curves in the outer disk. Most RCs of the other half are flat, and the small remainder are rising. Only 2 of the 41 galaxies are major mergers. Ten galaxies are demonstrably interacting with a small satellite, or a group, within about 1''–3'' from the central galaxy, and an additional five galaxies have projected small companions that might be satellites. The tidal forces of these small satellites in principle can affect the gas dynamics of the central galaxy in regions closest to the satellite, but we estimate that these perturbations are likely at the 10% level in most of these cases. Furthermore, in none of the cases with such nearby satellites is the overall RC asymmetric. Such perturbations and buckling/bending instabilities can drive warps of the main outer disks. Based on residual velocity maps (Paper II), deviations from circular rotation in a flat disk do occur in some of our RC41 sample but are not very large in amplitude, consistent with the fact that high-z disks are kinematically hot and geometrically thick, which in turn is expected to suppress or dampen the buckling and bending instability modes.

Tables D1 and D2 list the best-fit properties derived from our dynamical models of the RC41 sample. Columns 8 and 9 of Table D2 list the dark matter fractions at R_e obtained from the constrained posterior fitting (Appendix A.4) of the entire RC (out to 1.5–3.8 R_e), with 1σ uncertainties derived from the Markov Chain Monte Carlo (MCMC) analysis, which is discussed in more detail in Paper II. Figure 6 shows dark matter fractions, f_DM(R_e), as a function of disk circular velocity v_c(R_e) for z = 0 disks (upper-left panel), and passive, early-type galaxies from ATLAS-3D (top right), compared to our RC41 sample in two redshift slices (bottom panels) and to the results of the Illustris-TNG simulation (Lovell et al. 2018).

Zoom In Zoom Out Reset image size

Inverse Correlation between Galaxy Mass and Dark Matter Content. As is known from studies of low-redshift, star-forming disks, the dark matter fraction at the disk effective radius broadly decreases with increasing stellar, baryonic or halo mass, or circular velocity v_c(R_e) (Courteau et al. 2014 and references therein; Courteau & Dutton 2015). The gray line fit in the upper-left panel of Figure 6 (f_DM = 1–0.279 × (v_c −50.3 km s⁻¹)/100 km s⁻¹) shows this well-known inverse relationship. Dwarf galaxies are plausibly dark matter dominated, while the most massive, bulged disks near the Schechter mass have lower dark matter content and approach "maximal disks" (f_DM(max) < 0.28). The 15 SFGs in the lower redshift slice (z = 0.65–1.2) of our RC41 sample, as well as the average of the Wuyts et al. (2016) KMOS^3D sample (106 SFGs) in that redshift bin, also mostly fall on this inverse z ∼ 0 SFG correlation (lower-left panel of Figure 6).

In contrast, very few RC41 SFGs in the higher redshift slice (z = 1.2–2.45, lower-right panel of Figure 6) are found near the gray line. Seventeen of the 26 RC41 galaxies and the Wuyts et al. (2016) average in this redshift bin (92 SFGs) fall across a wide swath of low f_DM ≤ f_DM(max), with no or little dependence on v_c. This finding confirms and strengthens the earlier results of Wuyts et al. (2016), Genzel et al. (2017), and Lang et al. (2017), now with an almost seven times larger galaxy sample.

In the z = 1.2–2.45 redshift slice, more than 65% of the RC41 disks are baryon dominated within 1–3 R_e, and have dark matter fractions comparable to or less than maximal disks. Eleven of 17 SFGs (65%) in this region have f_DM ≤ 0.15, consistent with little or no dark matter, given the typical uncertainties. In contrast, at z = 0.65–1.2, only 4 of 15 (27%) have low dark matter fractions. The work of Courteau & Dutton (2015), Barnabe et al. (2012), and Dutton et al. (2013) indicate that at z = 0, the fraction of baryon-dominated, massive, bulgy disks in the same v_c and mass range is still smaller (<10%). Together, these findings suggest that cosmic time plays a role in setting the central dark matter fractions, in addition to mass or v_c. If correct, this conclusion could explain the inconsistencies between some of the recent studies at z ∼ 0.5–3, which we discussed in Section 1.3. For instance, Tiley et al. (2019b), Girard et al. (2018), and Di Teodoro et al. (2016, 2018) focus on galaxies at z ∼ 0.6–1.6 and with log(M_∗/M_⊙) ∼ 9.5–10.7. They find flat or even rising rotation curves, in excellent agreement with RC41 in the same part of parameter space. As found in z ∼ 0 galaxy studies, the dark matter fraction at one to several R_e is correlated with the shape of the RCs. In RC41, flat or dropping RCs yield modest or low dark matter fractions (0 $\to$ 0.5), while rising RCs signal larger dark matter fractions (up to 0.8).

Genzel et al. (2017) already pointed out that the low dark matter content within 1–3 R_e in many of their galaxies means that extrapolation of the dark matter mass to the virial scale with an NFW distribution results in unphysically large baryon to dark matter fractions averaged out to R_virial: m_baryon > f_b,c ∼ 0.17, the cosmic baryon fraction. This suggests that the low dark matter content in the outer disks cannot be extrapolated to the entire halo (Genzel et al. 2017; Lang et al. 2017). We will return to this point in Section 4.1.

Comparison of High-z Massive SFGs with z = 0 ETGs. Low dark matter fractions across a wide swath of v_c are characteristic of the z ∼ 0 ATLAS-3D passive, early-type galaxies (ETGs; upper-right panel of Figure 6; Cappellari et al. 2013, 2012; see Genzel et al. 2017). The massive z ∼ 2 SFGs are near or above the Schechter mass, M_S,* ∼ 10^10.7–10^10.9 M_⊙ (e.g., Ilbert et al. 2010, 2013; Peng et al. 2010). They are likely to quench and transition to the ETG population soon after z ∼ 2, or merge with another galaxy to form an even more massive ETG (Conroy et al. 2008; Genel et al. 2008). At least a fraction of the ATLAS-3D ETGs could be descendants of the massive z ∼ 2 SFGs. The low dark matter fractions of massive z = 0 ETGs could then already be set in their star-forming ancestors at z ∼ 1.5–2.5. As a cautionary note, the typically 10%–20% dark matter fractions within R_e of the ATLAS-3D ETGs in the upper-right panel of Figure 6 do assume a mass-dependent, bottom-heavy IMF (∼Chabrier/Kroupa at v_c ∼ 150 km s⁻¹, Salpeter at ∼240 km s⁻¹, super-Salpeter at ∼300 km s⁻¹; Cappellari et al. 2012). For a Chabrier/Kroupa IMF, dark matter fractions would go up to 25%–35% for the ETGs in the relevant v_c range (black upward-pointing arrow in the upper panels of Figure 6). We note that T. Mendel et al. (2020, in preparation) have observed a sample of massive z ∼ 1.4 ETGs as part of the KMOS-Virial survey. They find that essentially all of these ETGs, have <30% dark matter fractions, in agreement with our findings (see Figure 8).

In summary, the RC41 sample confirms and substantially strengthens the findings of Genzel et al. (2017) and Lang et al. (2017) that at z ∼ 2, massive SFGs with low or even negligible dark matter fractions within 1–3 R_e are common. With the almost seven times larger sample presented here compared to Genzel et al. (2017), we can now state that more than two-thirds of the massive SFGs at z = 1.2–2.45 are baryon dominated. Given the selection bias of RC41 to large galaxies (Figure 1), this fraction is likely a conservative lower limit. Such low dark matter, large massive SFGs exist also at z = 0–1, but they appear to be much rarer there. Similarly, low dark matter, baryon-dominated, passive, early-type galaxies at z = 0–2 may be descendants of the baryon-dominated z ∼ 2 population.

We now explore more quantitatively which galaxy properties correlate most strongly with the dark matter content discussed in the last section. Figure 7 shows representative correlation plots. We have also carried out a median-split correlation analysis and a broader principal component analysis (PCA), both of which we discuss in Appendix B. The different panels of Figure 7 show f_DM(R_e) on the vertical axis and the respective galaxy property on the horizontal axis. RC41 data are split into z = 0.65–1.2 (blue circles, 15 SFGs) and z = 1.2–2.45 (red circles, 26 SFGs) redshift slices. The large crosses denote the median uncertainties in these parameters. Best linear fit slopes (and 1σ errors) and the Pearson r² correlation coefficient are given in red. The thick gray line denotes binned data averages. Figure 7 shows that dark matter fractions at 1–3 R_e correlate most strongly with baryonic surface density within R_e, Σ_baryon (top right), with baryonic specific angular momentum λ × (j_baryon/j_DM) (top middle), and with bulge mass M_bulge (bottom right). The correlation with baryonic mass (bottom left), v_c, M_*, M_DM (bottom center), and galaxy size R_e (top left) are statistically not very significant (r² < 0.5). Median-split and PCA analyses give similar results (Appendix B). The PCA analysis also shows that there are no other significant correlations.

Zoom In Zoom Out Reset image size

The RC41 sample is biased toward relatively large systems to ensure well-resolved kinematics with subarcsecond resolutions. We are thus missing smaller and denser systems (Figure 1). Given the strong correlation of baryon fraction (=1 – dark matter fraction) with baryon density in the upper-right panel of Figure 7, we expect to find many more low dark matter, baryon-dominated systems in the missing population, consistent with the observations of van Dokkum et al. (2015), Wuyts et al. (2016), and Price et al. (2016, 2020), and the simulations of Zolotov et al. (2015).

By including the Wuyts et al. (2016) sample (Appendix C), we can extend the results to a much larger sample that better captures the entire z = 0.6–2.6 MS star-forming population and includes smaller, lower mass, and higher surface density galaxies. As discussed in the Introduction, Wuyts et al. (2016) have inferred the baryon fractions of 240 SFGs from KMOS^3D by comparing the dynamical mass obtained from the rotation curves in the inner disks (to ≤R_e) to the sum of stellar and gas masses, thus yielding f_baryon(R_e) = M_baryon/M_dynamical = 1 – f_DM(R_e). Most of the Wuyts et al. (2016) cubes are not as deep as those of RC41, and thus they do not have the extra constraint of the RC shape in the outer disk. The dark matter inference thus relies completely on the baryonic mass estimates based on stellar population synthesis modeling of the UV to near-IR SEDs (including the IMF; Wuyts et al. 2011a) and the gas scaling relations (Tacconi et al. 2018). Yet the Wuyts et al. (2016) sample provides a check of our results on a six times larger sample. In Figure 8, we show how the RC41 and Wuyts et al. (2016) results compare in the stellar mass–baryonic surface density plane (see upper-left panel of Figure 7). Keeping in mind the strong dependence of the Wuyts et al. (2016) dark matter values on systematic uncertainties of the input priors, the agreement of the two samples is remarkable. The Wuyts et al. (2016) sample indeed finds very low dark matter fractions at high surface density and many more low dark matter, baryon-dominated galaxies in the z = 1.2–2.6 slice than in the 0.6–1.2 slice. Both samples follow the same strong inversion correlation discussed in the last section. The massive ETGs at z = 1.4 (T. Mendel et al. 2020, in preparation) and z = 0 (Cappellari et al. 2012; Cappellari 2016) are located in the same region as the massive baryon-dominated z ∼ 2 SFGs (right panel of Figure 8).

Zoom In Zoom Out Reset image size

The continuous and dashed thick, light-blue curves in Figure 6 show the predictions of the Illustris TNG100 simulation at z = 0 (top panels) and z = 2, with DM profiles taken from a DM-only simulation and for a full hydro simulation, respectively (Lovell et al. 2018; M. R. Lovell 2020, private communication). Likewise, the gray curve in the right panel of Figure 8 shows the z = 2 TNG100 full hydro simulation with baryon–dark matter interactions in the Σ_baryon–f_DM plane. Broadly, the simulations capture the inverse correlation between f_DM and Σ_baryon (Figure 8) and provide a reasonable match to the observed dark matter fractions in z = 0 star-forming disks/spirals. The simulations also agree with the observational finding that dark matter fractions at a given v_c or Σ_baryon are greater at lower redshift than at high redshift. In more detail, there are significant differences between the simulations and observations. The baryon–dark matter interaction simulations over-predict dark matter fractions at z ∼ 2, especially at the high v_c, or high Σ_baryon tail. In Figure 6, only the pure dark matter simulations come close to the data. This is presumably because the baryon–dark matter interactions pull dark matter efficiently inward, which leads to higher dark matter concentrations at high z than in a dark matter only model. This is the effect of "adiabatic contraction" (AC; Blumenthal et al. 1986; Mo et al. 1998). The simulations could also miss essential physical processes, presumably on subgalactic scales, that keep the dark matter fractions as low as in dark matter only simulations. The "Magneticum" simulations (Teklu et al. 2018; K. Dolag et al. 2020, in preparation: thick brown line in the right panel of Figure 8) come much closer to the observed relation and sharp downturn of dark matter content with surface density. This is because this simulation includes strong AGN feedback.

Our most important finding of the last section is that 65% (RC41) or more (RC41 + sample of Wuyts et al. 2016) of massive star-forming disks at z ∼ 2 are baryon dominated within 1–3 R_e, and have dark matter fractions comparable to or less than maximal disks (〈f_DM〉 = 0.12 < f_DM(max) = 0.28). At z < 1.2, that fraction is 27%, and at z ∼ 0, less than 10% of the star-forming population appears to have such low dark matter fractions. Of these baryon-dominated, massive (〈M_baryon〉 ∼ 10¹¹ M_⊙) SFGs, 11 of 17 formally have f_DM ∼ $0\to$ 0.15, consistent with little or no dark matter, given the typical uncertainties. They are large (〈R_e〉 = 5.5 kpc) and have median baryonic angular momentum parameters, log(〈λ × (j_baryon/j_DM)〉) = −1.4, somewhat above the median of the overall population (−1.43, Burkert et al. 2016), as expected (Figure 1), but by much less than the rms scatter of the population (σ(log(λ × (j_baryon/j_DM)) ∼ 0.18). They are not "blue nuggets," proposed to be galaxies shortly after a rapid "compaction event," perhaps triggered by a major merger (Barro et al. 2013; Zolotov et al. 2015). While blue nuggets are more abundant at higher z, they are relatively rare and make up only about 10%–15% of the z ∼ 2 SFG population in the 10 < log(M_∗/M_⊙) < 11 mass range (Barro et al. 2013).

We have used the fitting functions of Moster et al. (2018, Appendix A3, Equation (A13); see also Behroozi et al. 2010, 2013; Moster et al. 2013; Burkert et al. 2016) as a separate estimate of the RC41 halo masses, given the redshift and the best-fit inferred RC41 stellar masses (column 3 of Table D1). The resulting dark matter masses are listed in column 4 of Table D1. We then computed ${m}_{d,{\rm{baryon}}}({\rm{Moster}})\,={M}_{{\rm{baryon}}}/{M}_{{\rm{DM}}}({M}_{\ast },z{)}_{{\rm{Moster}}}$ (column 8 of Table D1, filled brown triangles in Figure 9), which is the ratio of the observed baryonic mass in the disk (as obtained from the RC fitting) to the dark matter mass within R_virial, as inferred from the observed stellar mass and redshift and Equation (A13). We also computed ${m}_{d,{\rm{baryon}}}^{{\prime} }={M}_{{\rm{baryon}}}/{M}_{{\rm{DM}}}({\rm{NFW}},{f}_{{\rm{DM}}}({R}_{e}))$ (column 9 of Table D1, red and blue filled circles in Figure 9 for the usual two redshift slices), which is the ratio of the same observed baryonic mass in the disk, divided by the virial mass of the NFW halo (as obtained from the RC fitting). This halo has a dark matter fraction f_DM(R_e), which is identical with column 8 in Table D2. The thick horizontal gray line in Figure 9 shows the value that m_d would have if the cosmic baryon fraction, f_bc = 0.17, were all in the disk. Because in reality most of the baryons still reside in the halo, this line represents the upper physical limit of m_d. The triangles scatter around m_d,baryon ∼ 0.055, which is typical for the entire z ∼ 1–2.5 SFG population (Burkert et al. 2016). Triangles and circles agree with each other for f_DM(R_e) ≥ 0.35. For f_DM(R_e) < 0.3, however, the extrapolation to the virial radius of an NFW profile with the fitted f_DM(R_e) yields too high, or even completely unphysical, values for m_d'. This means that the assumption of an NFW dark matter profile on all radial scales cannot be correct in these cases.

Zoom In Zoom Out Reset image size

In our opinion, the most likely cause for the very low dark matter fractions in the majority of our z > 1.2 RC41 galaxies are deviations from the NFW profile at small radii. Flat or cored central dark matter distributions have been inferred for many local dwarf and spiral galaxies (e.g., Flores & Primack 1994; Moore 1994; Burkert 1995; McGaugh & de Blok 1998; de Blok & Bosma 2002; Marchesini et al. 2002; Kuzio de Naray et al. 2006; de Blok et al. 2008; Newman et al. 2012; Faerman et al. 2013; Adams et al. 2014; Oh et al. 2015), including the MW (Wegg et al. 2016; Portail et al. 2017).

To explore how deviations of the dark matter density distribution, ρ_DM, from the proposed original profile would affect dark matter fractions, we took Equation (A2), but now fitted the RC41 data with a combination of a bulge, a disk, and a modified NFW halo, with the additional free fit parameter α_inner(≥0), the central slope of the distribution. NFW has α_inner = 1, while a flat core (such as Burkert 1995) has α_inner = 0. We constrain the total mass for the modified NFW distribution by setting the integral of (A2) equal to the inverted Moster et al. (2018) estimate (column 4 in Table D1, Equation (A13)). Columns 10 and 11 of Table D2 give the fit values (and 1σ uncertainties) for α_inner.

Figure 10 shows the correlation between f_DM and these fitted values of α_inner. As expected, very low dark matter fractions are consistent with flat or cored dark matter distributions. Given the typical virial radii (column 3 in Table D2) and assumed concentration parameters (column 13 in Table D2), the transition to a cored distribution typically would occur at R_virial/c ∼ 24 kpc. Indeed, considering cored Burkert DM distributions (Burkert 1995) with variable α_inner and fitting for the core radius R_core as the free variable yields 〈 R_core〉 ∼ 25 kpc, consistent with the modified NFW analysis, albeit with large uncertainties. These core radii are also consistent with dark halo core scaling relations (Burkert 2015).

Zoom In Zoom Out Reset image size

We conclude that the low dark matter content of the massive high-z SFGs in the outer disks occurs naturally if the central cusp of the NFW profile is replaced by a cored or at least a less cuspy distribution than α_inner = −1. Replacing the Moster et al. (2018) scaling relations by those of Behroozi et al. (2013), or going from the concentration parameters listed in Table D2 to those one would derive from the Bullock et al. (2001b), Ludlow et al. (2014), or Dutton & Macciò (2014) fitting functions, does not qualitatively change these conclusions.

Several types of models might explain the presence of dark matter cores. The first explanation considers different dark matter particles, such as warm, fuzzy, or self-interacting dark matter (Hu et al. 2000; Spergel & Steinhardt 2000; Bode et al. 2001; Bertone et al. 2005; Calabrese & Spergel 2016; Hui et al. 2017; Pozo et al. 2020). The second favors a fundamental change of the law of gravity, MOND (Milgrom 1983). The third type of model considers interactions between baryons and dark matter in the inner halo and circumgalactic medium, where baryonic effects may act rapidly and efficiently (Governato et al. 2010; Macciò et al. 2012; Di Cintio et al. 2014; Chan et al. 2015; Peirani et al. 2017; Teklu et al. 2018; K. Dolag et al. 2020, in preparation).

On balance, a large abundance of baryon-dominated, dark matter cored galaxies at z ∼ 2, most strongly correlated with baryonic surface density, angular momentum, and central bulge mass, may be most naturally accounted for by the interaction of baryons and dark matter during the formation epoch of massive halos. Massive halos (log(M_halo/M_⊙) > 12) formed for the first time in large abundances in the redshift range z ∼ 1–3 (Press & Schechter 1974; Sheth & Tormen 1999; Mo & White 2002; Springel et al. 2005). At the same time, gas accretion rates were maximal (Tacconi et al. 2020). This resulted in high merger rates (Genel et al. 2008, 2009; Fakhouri & Ma 2009), very efficient baryonic angular momentum transport (Dekel et al. 2009; Zolotov et al. 2015), formation of globally unstable disks, and radial gas transport by dynamical friction (Noguchi 1999; Immeli et al. 2004; Genzel et al. 2008; Bournaud & Elmegreen 2009; Bournaud et al. 2014; Dekel & Burkert 2014). These processes enabled galaxy mass doubling on a timescale < 0.4 Gyr at z ∼ 2–3, and massive bulge formation by disk instabilities and compaction events on <1 Gyr timescales. However, central baryonic concentrations would naturally also increase central dark matter densities through adiabatic contraction (Barnes & White 1984; Blumenthal et al. 1986; Jesseit et al. 2002). For adiabatic contraction to be ineffective requires the combination of kinetic heating of the central dark matter cusp by dynamical friction from in-streaming baryonic clumps (El-Zant et al. 2001; Goerdt et al. 2010; Cole et al. 2011), with feedback from winds, supernovae, and AGNs driving baryons and dark matter out again (Dekel & Silk 1986; Pontzen & Governato 2012, 2014; Martizzi et al. 2013; Freundlich et al. 2020; K. Dolag et al. 2020, in preparation). Using idealized Monte Carlo simulations, El-Zant et al. (2001) demonstrated that dynamical friction acting on in-spiraling gas clumps can provide enough energy to heat up the central dark matter component and create a finite dark matter core (see also A. Burkert et al. 2020, in preparation). They argue that dark matter core formation in massive galaxies would require that clumps be compact, such that they avoid tidal and ram-pressure disruption, and have masses of >10⁸ M_⊙. Other idealized simulations (e.g., Tonini et al. 2006) confirm these results.

In full cosmological simulations, repeated and rapid oscillations of the central potential due to supernova and/or AGN feedback result in strong outflows and an expansion of the central collisionless component(s) (Read & Gilmore 2005; Mashchenko et al. 2008; Macciò et al. 2012; Pontzen & Governato 2012; Martizzi et al. 2013; Di Cintio et al. 2014; Chan et al. 2015; Read et al. 2016; Peirani et al. 2017; van der Vlugt & Costa 2019; K. Dolag et al. 2020, in preparation). However, these results are often based on galaxies at specific mass and redshift ranges and do not generally apply to the massive SFGs at z ∼ 0.5–2.5 discussed here. Nevertheless, the presence of giant star-forming clumps and the ubiquity of both stellar feedback-driven and AGN-driven outflows in the population of massive main-sequence galaxies at z ∼ 1–3 (e.g., Förster Schreiber et al. 2011, 2019) suggest that the above processes may be at work in RC41 galaxies.

As already mentioned, the results of the currently highest resolution, zoom-in cosmological hydro-simulations (such as Illustris-TNG50/100: Lovell et al. 2018; Nelson et al. 2019; Pillepich et al. 2019) qualitatively agree that high-z massive galaxies have lower dark matter fractions than those at z ∼ 0 but do not predict the "cored" dark matter distributions we infer from RC41 (Figure 6). Future very high spatial resolution cosmological simulations may properly capture dynamical friction and nuclear outflows on subcloud scales.

Assuming for simplicity that the growing virialized dark matter halo plus galaxy system starts with an initial state of a standard NFW dark matter distribution, the "final," cored state would require that dark matter from the cusp be partially removed to beyond the cusp's core radius. How large is this DM mass for the RC41 sample?

Here we compare again the difference in dark matter mass between the initial pure NFW halo (as estimated from Moster et al. 2018, Equation (A13)) and the DM inside R_e for the best-fit f_DM NFW models from Section 3. We then plot this "dark matter" deficit within R_e as a function of bulge mass, because bulge mass is one of the predictors of whether such a deficit (low dark matter fraction) occurs.

The results are listed in column 6 of Table D2 and shown in the left panel of Figure 11. Typical errors (large black cross) are large because we are computing differences of two uncertain numbers. Bulge masses are also more uncertain than baryonic masses. Nevertheless, within these uncertainties, the low dark matter fractions of about two-thirds of the RC41 galaxies are consistent with the hypothesis that on average dark matter equivalent to 30 (±10)% of the final baryonic bulge mass was removed from the galaxy core during bulge formation.

Zoom In Zoom Out Reset image size

Another approach is to exploit the fact that central dark matter column densities in NFW models at z ∼ 0.5–2.5 are fairly constant, 1–3 × 10⁸ M_⊙ kpc⁻² (column 4 of Table D2) and are only weakly dependent on z, λ, and M_∗. In the right panel of Figure 11, we use the N = 270 z = 0.6–2.6 star-forming disks with bulge masses inferred from multiband HST photometry (Lang et al. 2014) from the "full" KMOS^3D/SINS and zC-SINF sample (Wisnioski et al. 2015, 2019; Burkert et al. 2016) and plot them against the average dark matter surface densities within R_e (as estimated from the HST data), obtained from the Moster et al. (2018) scalings (colored distribution). The observed 〈Σ_DM(<R_e)〉−M_bulge data from the RC41 sample are marked as crossed green circles (column 4 of Table D2), clearly demonstrating the deficit compared to the Moster/Behroozi expectations. The filled brown triangles represent 〈Σ_DM(<R_e)〉 + 0.3 × M_bulge/πR_e² as a function of M_bulge. The distribution of the brown triangles is plausibly the same as that of the "full" sample or for NFW models.

Further support for a causal connection between the dissipative generation of dark matter cores during the first, gas-rich phase of massive galaxy formation comes from the HST morphologies of the RC41 SFGs. All five z > 2, dark-matter-rich or -dominant SFGs with f_DM(R_e) = 0.4–0.9 are clumpy rings without or with only small central bulges, while the baryon-dominated systems have very large bulges, with $\langle B/T\rangle \sim 0.5$ (column 17 of Table D1). At z < 1, only two of seven dark-matter-rich systems are ring dominated, while the remaining five have substantial bulges.

In Section 4.1, we noted that the RC41 SFGs are not "blue nuggets." This is strictly true in the sense that the blue nuggets when formed are compact, while most of the RC41 galaxies have extended disks/rings.  However, simulations suggest that after a compaction event a new, long-lived extended gas disk can re-form from inflowing high-angular-momentum gas and cold streams, whose inward mass transport is weakened by the presence of the massive central body (Dekel et al. 2020). Some of the RC41 galaxies, especially at z < 1.5, may be such galaxies.

Model Assumptions. As in our earlier work (Genzel et al. 2006, 2017; Burkert et al. 2016; Wuyts et al. 2016; Lang et al. 2017; Übler et al. 2017, 2018), we use forward modeling from a parameterized, input mass distribution. This mass distribution is assumed to be combination of a central unresolved bulge M_bulge (not radiating in the line emission tracer), a baryonic (gas plus stars) rotating disk M_disk(n_S, R_e, σ₀,i, pa) with Sérsic index n_S and disk effective radius R_e,¹⁵ inclined by an angle i relative to the plane of the sky, and with its major kinematic axis at position angle pa projected on the sky (pa positive east of north). The disk is assumed to have an isotropic, spatially constant velocity dispersion σ₀ (Wisnioski et al. 2015; Übler et al. 2019). For n_S = 1, exponential disks characteristic for most of RC41, the RC is given by Equation (A3). We assume a flat disk, with a z scale height determined by hydrostatic equilibrium (Spitzer 1942). We also assume that the disk effective radius R_e, which contains half the baryonic mass (gas and stars), is also the disk half-light radius (of the line emission tracer, Hα or CO).¹⁶ The third component is an extended dark matter halo. The circular velocity of the composite system at radius R then is

$$\begin{eqnarray}&&{v}_{\mathrm{circ}}^{2}(R)\,=\,{v}_{\mathrm{bulge}}^{2}(R)+{v}_{\mathrm{disk}}^{2}(R)+{v}_{\mathrm{DM}}^{2}(R).\end{eqnarray} \tag{ A1 }$$

For this model, we compute from Equation (A1) the three-dimensional circular velocity (v_c(x, y, z)) and line intensity (I(x,y, z)) distributions, rotate the model given input inclination and position angle of the line of nodes, and compute a 3D "model cube" I(x, y, vz) after integration of the rotated model distribution along the line of sight. We next convolve with the instrumental response function PSF(x, y, v_z) as appropriate for our measurements.

"Position–Velocity Cuts." We discuss the analysis of the RC41 sample in more detail in Paper II of this series (S. Price et al. 2020, in preparation). As already pointed out in the Introduction, here we focus our analysis on one-dimensional, software slit extractions (slit width δpp) of the beam-smeared, projected velocity distribution along the kinematic major axis (position angle pa). For these "position–velocity cuts" $I(p,{v}_{z}){| }_{{pa},\delta {pp}}$ we then compute from Gaussian fits (in apertures of length δp ∼ 015–03 along the slit) the velocity centroid, velocity dispersion, and line intensity as a function of p, for both the model and measurement cubes. We either use a constant software slit width (typically ∼1–1.5 FWHM of the data set) or a fanned slit (with 5°–10° opening angle), especially for low-inclination targets where the expected isovelocity contours fan out (e.g., van der Kruit & Allen 1978). In cases where OH sky emission lines affect the position–velocity cuts modestly, we clip the affected channels and interpolate. In cases of very strong OH sky-line contamination, we eliminate the galaxy from the sample. The 1D method is suitable for all RC41 galaxies, because most of the RC information is contained along this major axis (see Genzel et al. 2006, 2017).

Dark Matter Model. For the dark matter halo, we assume a Navarro et al. (1996, 1997, henceforth NFW) distribution of total mass M_DM (within R₂₀₀; see Mo et al. 1998) and total baryon to dark matter ratio m_db = (M_bulge + M_disk)/M_DM. The cold dark matter paradigm (Blumenthal et al. 1984; Davis et al. 1985; Blumenthal et al. 1986) predicts that the dark matter density distribution is described by a universal two parameter analytic function (NFW), alternatively (Einasto 1965). The dark matter density distribution in the NFW model is given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{DM}}(R)\,=\,\displaystyle \frac{{\rho }_{0}}{{\left(\tfrac{R}{{R}_{s}}\right)}^{{\alpha }_{\mathrm{inner}}}\,\times \,{\left[1+\left(\tfrac{R}{{R}_{s}}\right)\right]}^{3-{\alpha }_{\mathrm{inner}}}},\end{eqnarray} \tag{ A2 }$$

where R_s = R_virial/c is the scale radius of the distribution, and R_virial(z) is the virial radius of the DM distribution, within which the mean mass density is ∼200 times the critical density of the universe, and c is the halo's concentration parameter.

Because for dark matter only distributions the dependence of concentration parameter on halo mass is shallow (logarithmic exponent −0.05 to −0.12, Bullock et al. 2001b; Ludlow et al. 2014; Dutton & Macciò 2014) and the range of dark matter halo masses in RC41 is modest, we adopt for simplicity c ∼ 10.9 × (1 + z)^−0.83 as fixed inputs for our fits (an average of the published work, column 13 in Table D2), without a mass dependence. The scatter between these different theoretical descriptions is about ±0.08 dex.

The NFW distribution is cuspy with an inner logarithmic slope of α_inner = 1 (equivalent to α_Einasto = 1.8). At the scale radius R_s, the density distribution varies as R⁻² and in the outer halo as R⁻³. The RC of an NFW distribution has a positive slope from the center to R_s, which means that dark matter in the inner core is colder than the gas outside R_s. The radial distribution of the dark-matter-only circular velocity as a function of R is given by

$$\begin{eqnarray}&&\begin{array}{l}{v}_{\mathrm{DM}}^{2}={v}_{\mathrm{virial}}^{2}\times \left(\displaystyle \frac{{R}_{\mathrm{virial}}}{R}\right)\\ \,\times \,\displaystyle \frac{\mathrm{ln}\left(1+R/{R}_{s}\right)-\left(R/{R}_{s}\right)\left(1+R/{R}_{s}\right)}{\mathrm{ln}\left(1+c\right)-c/\left(1+c\right)},\end{array}\end{eqnarray} \tag{ A3 }$$

where v_virial is the DM circular velocity at R_virial. The circular velocity reaches a maximum value at 2.2 R_s, remains flat(ish) to about 0.4–0.5 R_virial, and falls off beyond. Integration of Equation (A2) from R = 0 to R_virial then yields the total dark matter mass M_DM.

For an infinitesimally thin, exponential disk (n_S = 1; Freeman 1970), the disk rotation and dark matter circular velocities are given by

$$\begin{eqnarray}&&\begin{array}{l}{v}_{\mathrm{disk}}^{2}(y=R/2{R}_{d})=4\pi G\times {{\rm{\Sigma }}}_{0}{R}_{d}\times {y}^{2}\\ \,\times \,\left[{I}_{0}\left(y\right){K}_{0}\left(y\right)-{I}_{1}(y){K}_{1}(y)\right],\end{array}\end{eqnarray} \tag{ A4 }$$

where I_n and K_n denote modified Bessel functions of order n.

We assume that the disk is pervaded by a (constant) velocity dispersion σ₀ (Burkert et al. 2010, 2016), such that the ratio of rotation velocity to dispersion is v_disk/σ₀. If the velocity dispersion is isotropic (Wisnioski et al. 2015), the disk then is geometrically thick (q = h_z/R_d ∼ σ₀/v_disk). For this thick disk case, we instead use the rotation curves derived for flattened spheroids of intrinsic axis ratio q = c/a as given by Equation (10) of Noordermeer (2008), solved numerically for a given q and n_S. For v_disk/σ₀ ∼ 4 typical at z ∼ 2, the disk rotation velocity at ∼R_1/2 is about 10% smaller than for a Freeman thin disk (Noordermeer 2008). Further, the presence of a pervasive velocity dispersion in an exponential disk creates an outward "pressure" gradient that decreases the centripetal force. Due to this "asymmetric drift," the overall rotation velocity becomes (for assumed isotropic and constant velocity dispersion σ₀)

$$\begin{eqnarray}&&{v}_{\mathrm{rot}}^{2}\,=\,{v}_{\mathrm{circ}}^{2}+2{\sigma }_{0}^{2}\times \displaystyle \frac{d\mathrm{ln}{\rm{\Sigma }}}{d\mathrm{ln}R}\,=\,{v}_{\mathrm{circ}}^{2}-2{\sigma }_{0}^{2}\times \left(\displaystyle \frac{R}{{R}_{d}}\right).\end{eqnarray} \tag{ A5 }$$

The impact of the asymmetric drift correction is two-fold and can be dramatic.

At R > R_e = 1.68 R_d (n_S = 1), the dark matter mass fraction is f_DM(R) = v_DM²(R)/v_circ²(R), and we can rewrite Equation (A1), for simplicity for B/T = M_bulge/(M_bulge + M_disk) ∼ 0, as

$$\begin{eqnarray}&&\begin{array}{l}{v}_{\mathrm{rot}}^{2}(R)\mathop{=}\limits^{{v}_{\mathrm{bulge}}=0}{v}_{\mathrm{disk}}^{2}\\ \,\times \left[\displaystyle \frac{1}{1-{f}_{\mathrm{DM}}}-3.36\times {\left(\displaystyle \frac{{\sigma }_{0}}{{v}_{\mathrm{disk}}}\right)}^{2}\times \left(\tfrac{R}{{R}_{e}}\right)\right].\end{array}\end{eqnarray} \tag{ A6 }$$

If the disk is not truncated, and if f_DM is < 1, there is formally a finite radius, at which the rotation velocity becomes small and mathematically even drops to 0 (at the same time the z scale height becomes very large), the system loses the character of a rotation-dominated disk and instead becomes a dispersion-dominated spherical system. For ${f}_{\mathrm{DM}}\to$0, this critical radius is ${R}_{\mathrm{crit}}/{R}_{e}\sim 4.8\times {\left({v}_{\mathrm{disk}}/4{\sigma }_{0}\right)}^{2}$. This drop is much faster than for a baryon-dominated but thin Freeman disk (e.g., Figure 2 of Genzel et al. 2017). At or near R_e, the asymmetric drift correction, that is, the second term within the brackets on the right side of Equation (A6), ∼−0.2. This means that if the RC is dropping significantly near and above R_e, dark matter fractions have to be small, f_DM(R_e) < 0.2.

To showcase the changes expected ab initio for RCs as a function of redshift, we show in Figure A1 two cases with the same baryonic mass (M_baryon = 8 × 10¹⁰M_⊙), same m_db = 0.05, same λ × (j_d/j_DM) ∼ 0.03, and same B/T = 0.3. The left galaxy at z = 0 is the MW with c ∼ 13, the right galaxy is D3a_15504 at z = 2.38 with c = 4. The MW model has negligible pressure effects, while D3a_15504 has σ₀ > 50 km s⁻¹ and v_c/σ₀ ∼ 4. For the MW case at z = 0, the RC is flattish without strong recognizable features of disk and halo. This insensitivity of the RC near R_e to the galaxy parameters is called the "disk–halo degeneracy" (or "conspiracy": van Albada et al. 1985).

Zoom In Zoom Out Reset image size

At high z, three things happen. First, for constant λ × (j_d/j_DM) R_e ∝ (1 + z)⁻¹ in ΛCDM, the baryonic disk is more compact and the peak disk plus rotation velocity correspondingly greater, for constant baryonic mass. Second, because of the smaller concentration of the halo (without considering AC), the halo moves outward and the combined overall rotation curve has a recognizable "baryonic peak" with outer drop. Third, this effect is strongly exacerbated by the pressure term, such that the inferred (and observed) rotation curve drops sharply, indicating the baryonic dominance. The "disk–halo degeneracy" thus can be more easily broken at high z.

Mo et al. (1998) assume that once a self-gravitating baryonic disk forms, the dark halo contracts adiabatically when the gravitational interactions between baryons and dark matter become important. If so, this counteracts the lower halo concentrations and higher baryonic densities of the ab initio model discussed above. However, feedback from supernovae, massive stars, and AGNs act to expand the halo. Burkert et al. (2010) argue from an analysis of the kinematics of high-z disks that the dark halos did not contract substantially during gas infall and disk formation. We consider both adiabatically contracted (AC) and non-contracted halos in our analysis, and average the resulting dark matter fractions (in columns 8 and 9 of Table D2). The inferred central dark matter fractions are not sensitive to AC for models with both free baryonic and free dark matter masses. For a given best-fit model, however, the effect of AC on the inferred central dark matter fraction is typically ∼+10%.

Tidal torque theory (Peebles 1969; White 1984) suggests that near the virial radius, the centrifugal support of baryons and dark matter is small and given by the angular momentum parameter,

$$\begin{eqnarray}&&\lambda \,=\,\displaystyle \frac{{\omega }_{\mathrm{virial}}}{{\omega }_{\mathrm{virial},\mathrm{cs}}}\,=\,\varepsilon \displaystyle \frac{{J}_{\mathrm{DM}}/{M}_{\mathrm{DM}}}{{R}_{\mathrm{virial}}\times \,{v}_{\mathrm{virial}}}\,\sim \,\displaystyle \frac{{J}_{\mathrm{DM}}\times {E}_{\mathrm{DM}}^{1/2}}{{{GM}}_{\mathrm{DM}}^{5/2}},\end{eqnarray} \tag{ A7 }$$

where ω = v_rot/R is the angular speed (v_rot is the rotational/tangential velocity) at R, and "virial" and "cs" stand for "at the virial radius" and "centrifugal support" (ω_rot,cs = (GM/R³)^1/2). The constant is $\sim \sqrt{2}$, J and j are the total and specific (j = J/M) angular momenta, and E ∼ GM²/R is the absolute value of the total gravitational energy. Building up on earlier work by Peebles (1969) and Barnes & Efstathiou (1987), simulations have shown that tidal torques generate a universal, near-log-normal distribution function of halo angular momentum parameters, with 〈λ〉 = 0.035–0.05 and a dispersion of ±0.2 in log (Bullock et al. 2001a; Hetznecker & Burkert 2006; Bett et al. 2007; Macciò et al. 2007).

If the baryons are dynamically cold or can cool rapidly after shock heating at R_virial, they are transported inwards and form a centrifugally supported disk of (exponential) radial scale length R_d, given by, e.g., Fall (1983) and Mo et al. (1998),

$$\begin{eqnarray}&&{R}_{d}=\displaystyle \frac{1}{\sqrt{2}}\left(\displaystyle \frac{{f}_{{Jd}}}{{m}_{d}}\right)\times \lambda \times {R}_{\mathrm{virial}}=\displaystyle \frac{1}{\sqrt{2}}\left(\displaystyle \frac{{j}_{d}}{{j}_{\mathrm{DM}}}\right)\times \lambda \times {R}_{\mathrm{virial}}.\end{eqnarray} \tag{ A8 }$$

Here, m_d = M_d/M_DM is the ratio of the baryonic disk mass to that of the dark matter halo and f_Jd is the fraction of the total dark halo angular momentum in the disk, J_d = f_Jd J_DM. In the classical literature, it has generally been assumed that the baryonic angular momentum is conserved between the virial and disk scale, such that j_d ∼ j_DM (e.g., Fall 1983; Mo et al. 1998). More recent simulations find that various competing processes between the outer halo and inner disk scales can lead to up or down variations of j_d/j_DM by factors of 3 (Übler et al. 2014; Danovich et al. 2015; Teklu et al. 2015).

Changes of Angular Momentum due to Turbulent Pressure and Deviations from Exponential Distributions. The observations provide an estimate of σ = σ₀ at ∼2–2.5 R_1/2, which we adopt as the characteristic dispersion everywhere in the disk (see last section). According to Equation (A6), this isothermal disk has a finite "truncation" radius R_max/R_d = 0.5 × (v_circ/σ₀)² ∼ 2–15, where v_rot = 0. For kinematically cold disks with large ratios of rotation to dispersion, v_rot = v_circ and the solution approaches the constant value R_1/2 = 1.68 × R_d. Disks with larger velocity dispersions can, however, be strongly dispersion truncated with half-mass radii that can become even smaller than R_d. A convenient approximation is (Burkert et al. 2016)

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{1/2}}{{R}_{d}}=1.7{\left(\displaystyle \frac{{v}_{1/2}}{{\sigma }_{0}}\right)}^{2.7}\times {\left(1+{\left(\displaystyle \frac{{v}_{1/2}}{{\sigma }_{0}}\right)}^{2.7}\right)}^{-1}.\end{eqnarray} \tag{ A9 }$$

Here, v_1/2 = v_rot(R_1/2).

Appendix B of Burkert et al. (2016) reformulates Equations (A7) and (A8) in the framework of the Mo et al. (1998) NFW model to yield the following fitting function for the baryonic angular momentum parameter λ × (j_d/j_DM):

$$\begin{eqnarray}&&\begin{array}{l}\lambda =f\times 5.16\times {10}^{-5}\times \left(\displaystyle \frac{{j}_{\mathrm{DM}}}{{j}_{d}}\right)\times \left(0.0153c+1\right)\\ \times \,\left(\displaystyle \frac{{R}_{d}\times {v}_{1/2}}{\,\mathrm{kpc}\times \,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)\times {\left(\displaystyle \frac{H(z)}{100\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{kpc}}^{-1}}\right)}^{1/3}\\ \times \,{\left(\displaystyle \frac{{M}_{\mathrm{baryon}}/{m}_{d}}{{10}^{12}{M}_{\odot }}\right)}^{-2/3}\\ f=1\,\,\mathrm{for}\,\left({v}_{1/2}/{\sigma }_{0}\right)\gt 10\\ f=0.5\left[1+\cos \left(\pi \left({\sigma }_{0}/{v}_{1/2}-0.1\right)\right)\right]\\ \mathrm{for}\,\,10\geqslant \left({v}_{1/2}/{\sigma }_{0}\right)\geqslant 1.4\\ f=0.35\,\,\mathrm{for}\,\left({v}_{1/2}/{\sigma }_{0}\right)\lt 1.4.\end{array}\end{eqnarray} \tag{ A10 }$$

Here, c is the halo concentration parameter, c = R_c/R_virial (Equation (A2), column 13 of Table D2).

So far, we have assumed that the surface density distribution of baryons is exponential (n_S = 1). This is motivated by the CANDELS 3D-HST analysis of the H-band light of the reference population of the SFGs (Wuyts et al. 2011b; Bell et al. 2012; Lang et al. 2014). Twenty-five of the RC41 SFGs have n_S = 0.9–1.1. However, five of the RC41 galaxies have n_S = 1.2–1.3 and 10 have n_S = 0.2–0.6. Romanowsky & Fall (2012) and Burkert et al. (2016) analyzed the combined impact of turbulence (i.e., asymmetric drift) and deviations from exponential distributions on the specific angular momentum of baryonic disks. They cast the results in the form of a fitting function depending on the parameters $x=\mathrm{log}{n}_{S}$ and $y\,=\mathrm{log}({\sigma }_{0}/{v}_{1/2})$, such that

$$\begin{eqnarray}&&\begin{array}{l}k({n}_{S},y)={j}_{d}({n}_{S},y)/\left({v}_{1/2}\times {R}_{1/2}\right),\,\,\mathrm{with}\\ {\rm{\Sigma }}(R,{n}_{S})=\exp \left(-{b}_{{n}_{S}}{\left(-\displaystyle \frac{R}{{R}_{1/2}}\right)}^{1/{n}_{S}}\right)\,\mathrm{and}\\ {b}_{{n}_{S}}=2{n}_{S}-1/3+0.009876/{n}_{S},\,\,\mathrm{and}\\ \mathrm{log}k(x,y)=-0.082+0.091\\ \times \,x-(0.06+0.244\times y)\times {x}^{2}-0.168\times y.\end{array}\end{eqnarray} \tag{ A11 }$$

Equation (A11) fits all combined data in the interval x = −0.7 to 0.7 and y = −1.2 to −0.15 to better than ±0.03 dex. For the relevant range in x and y, the inferred values of k(x, y) vary from ∼1 to ∼1.75, where for a thin exponential disk k(0, −$\infty$) = 1.19. These corrections tend to slightly decrease λ × (j_d/j_DM) for SFGs at the low mass tail, and slightly increase λ × (j_d/j_DM) for SFGs at the high mass end of our sample. These corrected baryonic angular momentum parameters of the disks are listed in Column 7 of Table D2 and used for the top central panel in Figure 7.

Several groups have inferred mean stellar mass to dark matter halo mass ratios as a function of cosmic time and mass from ranked matching of observed stellar masses of galaxies with computed dark matter masses of the same abundance (both as a function of z and mass; see Behroozi et al. 2010, 2013; Moster et al. 2013, 2018). These analyses give the ratio of M_*/M_DM, or M_∗, as a function of z and M_DM. One can use their fitting functions to in turn compute M_DM if z and M_∗ are given, taking into account the conditional probability distributions. That is, $P({M}_{{\rm{DM}}}| {M}_{\ast })\,=P({M}_{\ast }| {M}_{{\rm{DM}}})\ast P({M}_{{\rm{DM}}})/P({M}_{\ast })$, where $P({M}_{\ast })\,=\int d{M}_{{\rm{DM}}}\,\times P({M}_{\ast }| {M}_{{\rm{DM}}})$.

Here we use stellar masses either from SED fitting or inferred from the best-fit baryonic mass together with a fixed gas fraction. For Moster et al. (2018), we have data for individual galaxy–halo pairs at all redshifts, and we can directly fit the halo mass–galaxy mass relation. We then find (B. Moster 2020, private communication):

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}{M}_{\mathrm{DM}}({M}_{\odot })\,=\,-0.301+\mathrm{log}{M}_{* }+a+\mathrm{log}(\delta {M}^{b}+\delta {M}^{c}),\\ \mathrm{with}\\ \delta {\rm{M}}={M}_{* }({M}_{\odot })/3.981\times {10}^{10}({M}_{\odot })\\ a=1.507-0.124\times (z/(1+z)),\\ b=-0.621-0.059\times (z/(1+z)),\\ c\,=\,1.055+0.838\times (z/(1+z))-3.083\times {\left(z/(1+z)\right)}^{2}.\end{array}\end{eqnarray} \tag{ A12 }$$

which works well for z > 0.5. Using the corresponding relations for Behroozi et al. (2013) gives very similar results. In the RC41 region (log M_* ∼ 10–11.3 and z = 0.65–2.4), we find $\langle \mathrm{log}{M}_{{\rm{DM}}}({\rm{M}}18)-\mathrm{log}{M}_{{\rm{DM}}}({\rm{B}}13) \rangle {| }_{z,\mathrm{log}{M}_{\ast }}=-0.02\ldots -0.1$. With M_DM in hand, it is then possible to compute the dark matter mass distribution M_DM(R) and surface density distribution Σ_DM(R) = M_DM(R)/πR² with Equations (A1) and (A3) using an NFW or modified NFW distribution (i.e., with α_inner < 1).

Marquardt–Levenberg gradient χ² fitting. We carried out two independent analyses of the RC41 data set. First, we used a classical Marquardt–Levenberg gradient χ² minimization of the p−v(p)–δv(p) and p–σ(p)–δσ(p) data points at offset p from the center, along the kinematic major axis, and compared to the extraction of our 3D models at the same positions and smeared to the same angular and velocity resolution. We also used as additional qualitative constraints the p–I(p)–δI(p) measurements of the Hα or CO line but did not include them in most but a few galaxies in the formal fitting because the Hα-light and mass distributions can be quite different in the center if a low-star-forming bulge is present. As the data are sampled every δp = 015–03, we obtain, depending on angular resolution and whether or not adaptive optics was used, between 1.5 and 3.5 samples per FWHM angular resolution. The total number of data points from the two velocity and velocity dispersion cuts per galaxy then varies between 22 and 60 so that we effectively have 20–30 independent data points (in a Nyquist-sampled sense). To characterize the mass distribution and kinematics, the minimum model parameters are M_baryon(disk) = M_bulge + (M_gas + M_*)_disk, B/T = M_bulge/(M_gas + M_*)_disk, σ₀ and M_{DM (NFW)}, for a total of four free fit parameters. In this case R_e, n_S (disk), i (inclination)_disk, q_disk, c_halo, adiabatic contraction AC_halo, are all input parameters constrained/given by ancillary information (mainly HST imaging data, or DM simulations). In some cases, it is necessary to include R_e, and in some rare cases also i (typically for face-on SFGs) as formal fit parameters, thus increasing the number of fit parameters to five to six. This leaves 15–25 degrees of freedom. In addition, multiband photometry analysis gives a prior on M_* (e.g., Wuyts et al. 2011a), and scaling relations give a prior on M_gas (e.g., Scoville et al. 2017; Tacconi et al. 2018).

MCMC Fitting. Our second approach used a full Bayesian MCMC investigation of the physical parameters. For the present paper, we focus on obtaining uncertainties for the key physical parameters, the dark matter fraction and the total baryonic mass. While the majority of the uncertainty on f_DM is captured when fitting only to f_DM and log M_baryon to mirror the χ² analysis, we perform MCMC fits with f_DM(R_e), log M_baryon, B/T, σ₀, and if needed, R_e free. The other parameters are fixed to the adopted ancillary parameter values (or the adjusted inclination value). We use flat priors on f_DM(R_e) and σ₀, and Gaussian priors on log M_baryon, B/T, and R_e using prior information from ancillary information (but using kinematically derived prior centers for B/T and R_e, as the light and mass distributions in galaxies can differ). Paper II discusses the methodology, results, and uncertainties in finding the best fitting values for the basic fit parameters and for the MCMC analysis.

The two methods broadly agree quite well, $\langle {f}_{{\rm{DM}}}({\chi }^{2}) \rangle \,- \langle {f}_{{\rm{DM}}}({\rm{MCMC}}) \rangle =-0.11(\pm 0.18)$, where the number in the bracket is the median 1σ rms uncertainty from the MCMC analysis. The MCMC analysis does give (marginally) higher DM fractions, mainly for galaxies of low dark matter content in the χ²-technique. The most likely explanation for the small systematic trend of −0.1 in f_DM is the treatment of priors in the two analyses. For the χ² analysis, we used a top-hat distribution of priors, while in the MCMC analysis, we used a Gaussian distribution for the prior on the baryonic mass. In the MCMC fitting, there is a trend to favor lower baryonic masses, a number of them more than one unit of dispersion below the prior maximum. This offset then increases the allowable dark matter fraction. The same effect is seen in Wuyts et al. (2016), where their Figures 5, 7, and 9 show a significant number of unphysical values of log (M_baryon/M_dyn) > 0 (and up to 0.5).

In the following, we use the values of f_DM obtained from the χ² fitting as our final estimators and adopt the MCMC uncertainties as the more robust estimators of uncertainties (Paper II).

Dependence of dark-matter fractions on parameter correlations. Figure A2 gives four examples of the dependence of total χ² on different galaxy parameters, for the NOEMA-CO data set in the PHIBSS1 galaxy EGS_13035123 (z = 1.12; Tacconi et al. 2010). These figures show that σ₀, R_e, and B/T are constrained to ±15%–20%, ±6%, and ±10%–15% (1σ) in this galaxy. The χ² dependence of f_DM(≤R_e) is shallower, and dark matter fractions can rarely be better determined than ±0.1 in absolute units or ±20%–40% in fraction. This fractional uncertainly is increased by at least another ±0.05 if it is unknown whether AC is effective (see discussion in Methods of Genzel et al. 2017).

Zoom In Zoom Out Reset image size

In terms of covariance, we find that increasing B/T leads to increased f_DM(R_e), which can be understood in the sense that a higher baryonic mass fraction in the central bulge decreases the relative contribution of the baryons to v_c at R_e. When changing R_e, the effects are less definite, but for the majority of cases, we find that increasing R_e again leads to increased f_DM(R_e). For these cases, this result can be understood in the sense that a larger R_e distributes the baryonic mass onto a larger disk (that is, less compact), leading to less relative contribution of the baryons to v_c at R_e. When increasing or decreasing the best-fit R_e or B/T by their uncertainties given in Figure A2, the resulting changes in f_DM are all below ∼0.13.

The inferred dark matter fraction depends on the stellar initial mass function (IMF) adopted. We used a Chabrier (2003) IMF for the input stellar priors. A bottom-heavy IMF, as proposed by van Dokkum & Conroy (2010) and Cappellari et al. (2012) for very massive passive galaxies, would imply a higher stellar mass for the prior information. A Salpeter IMF (compared to the Chabrier IMF) yields a 0.23 dex greater (M/L*) ratio than the Chabrier IMF. This change in the input prior may push the inferred dark matter fractions to lower values. However, at z ∼ 1–2.5, average dense gas to stellar mass ratios in the inner star-forming disks of SFGs are ∼0.5–1, thus lowering the impact of bottom-heavy IMFs correspondingly. However, our fitting method directly fits for the total baryonic mass, without regard to the form of the IMF. Thus, modulo changes to the priors, assuming a more bottom-heavy IMF would not change the dark matter fractions we find in this study.

The dark-matter fraction also depends on the mass distribution of the halo. If adiabatic contraction in the inner halo on 1–3 R_e is effective, dark matter fractions go up correspondingly, typically by 10% when using the Mo et al. (1998) AC recipes. Recall that our values for f_DM(R_e) in columns 8 and 9 of Table D2 are averages between AC on and off.

As mentioned above, we adopt a simple scaling of concentration parameter as a function of redshift only, which is an average of published work (Bullock et al. 2001b; Dutton & Macciò 2014; Ludlow et al. 2014). In reality, concentration parameters in simulations are mass dependent ($\mathrm{log}c\,\sim -0.05\ast \mathrm{log}{M}_{{\rm{DM}}}$) and show considerable scatter (e.g., Dutton & Macciò 2014; Ludlow et al. 2014). Lower concentration parameters (at higher z, higher mass, etc.) result in somewhat lower f_DM(R_e): in most cases, the baryonic mass and effective radius remain largely unchanged over a wide range of assumed concentration parameters, although reducing c (down to c = 2) can increase the inferred M_baryon by up to ∼0.2 dex and increase the effective radius by a few kiloparsecs. For instance, at z ∼ 2.3, an extreme change from c ∼ 4 to c ∼ 2 results in δf_DM = −0.05, which can be significant for high-mass, low dark matter systems. But generally, varying the concentration parameter does not have a large impact on the resulting fit parameters, including the dark matter fractions. For the galaxy displayed in Figure A2, large changes in c have only a negligible effect on f_DM(R_e) and χ², as shown in the upper-left panel (no AC is considered here). The inset shows the probability distribution function for halo concentration parameters based on the EMERGE model (Moster et al. 2018, 2019) for central halos at z = 1.1 with 10.6 < log(M_*/M_⊙) < 11.2 and 20 < SFR < 150 (M_⊙ yr⁻¹), where stellar mass and star formation rate include observational uncertainties.

In Paper II, we investigate the value of using the full 2D velocity and velocity dispersion maps, or, even more ambitiously, fitting the data cubes and full spaxel line profiles. Genzel et al. (2017, Methods section) have shown examples of the 2D fitting for a few well-resolved disks. Paper II shows that for most of the RC41 data sets most of the information is contained in the major axis cuts, and additional constraints (e.g., on inclination) from 2D data are only effectively available only for a few (mostly AO assisted) data sets. Spaxel profile fitting (or parameterizing deviations from Gaussianity with h3 and h4 (Naab & Burkert 2003; Emsellem et al. 2004) can be very helpful for constraining the inner disk/bulge kinematics and structure, and deviations from circular motions. For the outer disk kinematics, which is the focus of the present paper, such an additional analysis effort is not required, because the profiles are close to Gaussian shape in most galaxies.

Here we quantified trends between pairs of properties with the Pearson correlation coefficient and examined all properties simultaneously through PCA. We included as possible vectors the redshift, the masses characterizing each galaxy (the baryonic mass, the dark matter mass at the virial radius, and the central bulge mass), the effective disk size, and the disk scale angular momentum parameter of the baryons, the potential well depth at the effective radius (i.e., the circular velocity), the baryonic mass surface density within the effective radius, and the inner slope of the dark matter halo profile α_inner. We analyzed the RC41 sample alone and in combination with the Wuyts et al. (2016) "W16" sample (excluding from the latter 14 galaxies already in RC41). For the RC41+W16 case, the α_inner is omitted as it was derived only for the RC41 objects, and the W16 set is restricted to the subset of 133 galaxies that have a stellar bulge mass estimate (while all RC41 galaxies have total stellar+gas bulge masses determined from dynamical modeling).

Figure B1 summarizes the results of the PCA analysis. For the RC41 sample, 72% and 83% of the total cumulated fractional variance are captured by two and three principal components, respectively. For the four times larger but less well-constrained RC41+W16 sample (174 galaxies), 67% and 79% of the variance are captured by two and three principal components. In the 2D correlation matrices on the left, we sorted the galaxy parameters from top to bottom going from the strongest correlation to the strongest anticorrelation with f_DM(<R_e). Nearly identical results are found using the Spearman rank correlation coefficient instead of the Pearson coefficient. The loading plots from the PCA in the space of the first two principal components (PC1 and PC2), shown in the right-hand panels, reveal that PC1 is strongly coupled to Σ_baryon and various mass estimates, and PC2 is most tightly associated with galaxy size.

Zoom In Zoom Out Reset image size

Focusing on f_DM(<R_e), the correlation and PCA analyses clearly show that the DM fraction averaged over the disk within R_e is most strongly correlated with Σ_baryon and λ × (j_baryon/j_DM). For the RC41 sample, we also find a strong correlation with M_bulge, which becomes much weaker when including the W16 objects, perhaps because the bulge masses are based on stellar light measurements from HST data and are not sensitive to gas-rich or highly extincted bulges. The correlations with all other parameters are weak, including redshift, physical size, and v_c. A more significant trend is seen in RC41+W16 with redshift, presumably because this sample is larger and has many more SFGs in the lower redshift slice, z = 0.6–12, than RC41. Expanding the PCA and correlation analysis to all 226 W16 galaxies not in RC41 (for a total of 267) and excluding M_bulge has little impact on the results for all other properties, with correlation coefficients changing by 0.05 typically.

We also split the sample into the lower and upper half around the median of each parameter, and then evaluated how significant any difference was in terms of the combined uncertainties of the mean. We carried out this analysis for the RC41 and Wuyts et al. (2016) samples, for the same parameters as in Appendix B.1. Tables B1 and B2 summarize the results.

Table B1.  Median Dark Matter Fractions in Two Bins for a Given Parameter

Parameter	Value	St-dev	$\langle {f}_{{\rm{DM}}} \rangle$	Uncertainty	Significance
Σbaryon	8.47	0.37	0.42	0.03	⋯
	9.10	0.30	0.14	0.02	−7.5
λ * jd/jDM	−1.41	0.32	0.15	0.02	⋯
	−1.19	0.40	0.42	0.03	7.2
αinner	0.09	0.33	0.15	0.01	⋯
	1.10	0.40	0.42	0.03	7.1
Re	4.60	0.32	0.15	0.01	⋯
	7.10	0.44	0.34	0.03	5.3
Mbulge	10.20	0.37	0.34	0.03	⋯
	10.73	0.23	0.19	0.02	−3.9
Mbaryon	10.63	0.28	0.32	0.02	⋯
	11.12	0.18	0.20	0.02	−3.5
z	0.92	0.39	0.32	0.02	⋯
	2.20	0.36	0.22	0.02	−2.9
MDM	11.87	0.33	0.30	0.02	⋯
	12.30	0.22	0.22	0.02	−2.6
vc	194.000	0.321	0.250	0.022	⋯
	317	0.199	0.225	0.027	−0.7

Note. Columns 2 and 3 give the median parameter values (and their standard deviations) in the two bins, the top one being the half of the 41 RC41 galaxies below and the bottom one being the half above the median of a given parameter. Columns 4 and 5 give the median dark matter fractions (and their 1σ uncertainty) for the two bins. Column 6 gives the significance of the difference between the two bins (upper minus lower), in units of the combined uncertainty. The table is ordered top to bottom in decreasing level of significance.

Download table as:  ASCIITypeset image

Table B2.  Same as Table B1, but Now for the 240 Galaxies of the Wuyts et al. (2016) Sample

Parameter	Value	$\langle {f}_{{\rm{DM}}}({R}_{e}) \rangle$	Uncertainty	Significance
λ * jd/jDM	0.012	0.13	0.03	⋯
	0.030	0.59	0.03	11
Σbaryon	8.47	0.59	0.03	⋯
	9.06	0.24	0.03	−8.2
z	0.91	0.55	0.03	⋯
	2.19	0.26	0.03	−6.8
Mbaryon	10.40	0.52	0.03	⋯
	11.00	0.29	0.03	−5.4
MDM	11.69	0.49	0.03	⋯
	12.19	0.29	0.03	−4.7
Re	2.85	0.32	0.03	⋯
	5.04	0.51	0.03	4.5
Mbulge	9.36	0.46	0.03	⋯
	10.32	0.29	0.03	−4
vc	184	0.38	0.03	⋯
	292	0.50	0.03	2.8

Note. Note that only 139 SFGs of the Wuyts et al. (2016) sample have stellar bulge masses, estimated from HST imagery. In contrast, for the RC41 sample, we estimate the stellar and gas mass content in the bulge region. The Wuyts et al. (2016) bulges thus are lower limits and strongly effected by extinction.

Download table as:  ASCIITypeset image

The results in Tables B1 and B2 are in excellent agreement with those obtained in the correlation and PCA analyses in Appendix B.1.

To summarize, we have explored with three methods the significances of the correlation of the dark matter (or baryon) fraction averaged over the disks of z = 0.6–2.6 star-forming galaxies. All three agree that this quantity empirically is most strongly correlated with baryonic surface density, baryonic angular momentum parameter, and, somewhat less convincingly, with the total mass of the central bulge. While these findings may be valuable for understanding the origin of the low dark matter fractions in a surprisingly large number of high-mass SFGs at z ∼ 2, many of the parameters we have studied are interrelated, and it is difficult, and sometimes misleading, to extract causation from empirical correlations (see Lilly & Carollo 2016).