The MOSDEF Survey: Kinematic and Structural Evolution of Star-Forming Galaxies at $1.4\leq z\leq 3.8$

TL;DR

This study analyzes ionized gas kinematics of 681 high-redshift star-forming galaxies to understand their structural evolution, baryonic content, and dynamical properties, revealing trends with mass, star formation, and surface density.

Contribution

It provides the first large-scale kinematic analysis of galaxies at $z\,1.4$ to $3.8$, combining models accounting for misalignments with structural data to constrain galaxy dynamics and baryon fractions.

Findings

$V/\sigma_{V,0}$ increases with stellar mass and decreases with sSFR.

Galaxies are baryon dominated within their effective radii at high redshift.

Baryon fractions within $R_E$ tend to decrease over cosmic time, often exceeding 100%.

Abstract

We present ionized gas kinematics for 681 galaxies at $z\sim 1.4-3.8$ from the MOSFIRE Deep Evolution Field survey, measured using models which account for random galaxy-slit misalignments together with structural parameters derived from CANDELS Hubble Space Telescope (HST) imaging. Kinematics and sizes are used to derive dynamical masses. Baryonic masses are estimated from stellar masses and inferred gas masses from dust-corrected star formation rates (SFRs) and the Kennicutt-Schmidt relation. We measure resolved rotation for 105 galaxies. For the remaining 576 galaxies we use models based on HST imaging structural parameters together with integrated velocity dispersions and baryonic masses to statistically constrain the median ratio of intrinsic ordered to disordered motion, $V/\sigma_{V,0}$. We find that $V/\sigma_{V,0}$ increases with increasing stellar mass and decreasing specific SFR (sSFR). These trends may reflect marginal disk stability, where systems with higher gas fractions have thicker disks. For galaxies with detected rotation we assess trends between their kinematics and mass, sSFR, and baryon surface density ($\Sigma_{\mathrm{bar},e}$). Intrinsic dispersion correlates most with $\Sigma_{\mathrm{bar},e}$ and velocity correlates most with mass. By comparing dynamical and baryonic masses, we find that galaxies at $z\sim 1.4-3.8$ are baryon dominated within their effective radii ($R_E$), with Mdyn/Mbaryon increasing over time. The inferred baryon fractions within $R_E$, $f_{\mathrm{bar}}$, decrease over time, even at fixed mass, size, or surface density. At fixed redshift, $f_{\mathrm{bar}}$ does not appear to vary with stellar mass but increases with decreasing $R_E$ and increasing $\Sigma_{\mathrm{bar},e}$. For galaxies at $z\geq2$, the median inferred baryon fractions generally exceed 100%. We discuss possible explanations and future avenues to resolve this tension.

Tables

Table 1. Table 1: Kinematic variables

$V_{2.2}$; $V\left(R_E\right)$	Intrinsic rotation velocity at $r_{2.2}=2.2r_s\approx 1.3R_E$; $R_E$
${\sigma }_{V,0}$	Intrinsic galaxy velocity dispersion (assumed to be const.)
${\sigma }_{V,1D,\mathrm{obs}}$	Directly observed 1D integrated velocity dispersion
${\sigma }_{V,\mathrm{obs}}$	Integrated 1D velocity dispersion, corrected for inst. res.
${\sigma }_{V,\mathrm{model}}$	Kinematic model integrated 1D velocity dispersion

Table 2. Table 2: Median kinematic properties, Resolved/Aligned galaxies with detected rotation

			Median Properties
Redshift	Bin	$N$	$⟨z⟩$	$⟨{\log }_{10}\left(M_{\ast }\right)⟩$	$⟨{\log }_{10}\left(\mathrm{sSFR}\right)⟩$	$⟨{\log }_{10}\left(R_E\right)⟩$	$⟨{\left(V/{\sigma }_{V,0}\right)}_{R_E}⟩$	$⟨{\left(V/{\sigma }_{V,0}\right)}_{2.2}⟩$	$⟨\left|V_{2.2}\right|⟩$	$⟨{\sigma }_{V,0}⟩$
			—	$\left[{\log }_{10}\left(M_{\odot }\right)\right]$	$\left[{\log }_{10}\left({\mathrm{yr}}^{-1}\right)\right]$	$\left[{\log }_{10}\left(\mathrm{kpc}\right)\right]$	—	—	$\left[\mathrm{km}/s\right]$	$\left[\mathrm{km}/s\right]$
	$[8,9.75]$	5	1.59	9.49	-9.09	0.32	${1.9}_{-1.1}^{+0.8}$	${2.5}_{-1.5}^{+1.6}$	${122}_{-86}^{+1}$	${42}_{-12}^{+5}$
	${\log }_{10}\left(M_{\ast }\right)=[9.75,10.25]$	22	1.54	10.04	-8.90	0.54	${1.6}_{-0.3}^{+0.2}$	${2.0}_{-0.5}^{+0.2}$	${131}_{-29}^{+4}$	${67}_{-11}^{+2}$
	$[10.25,12]$	14	1.53	10.69	-9.16	0.66	${2.8}_{-0.7}^{+0.4}$	${3.4}_{-0.9}^{+0.5}$	${197}_{-16}^{+41}$	${70}_{-10}^{+9}$
	$[-11,-8.9]$	24	1.53	10.15	-9.13	0.56	${1.7}_{-0.3}^{+0.5}$	${2.2}_{-0.4}^{+0.6}$	${134}_{-14}^{+17}$	${67}_{-11}^{+2}$
$1.3\le z\le 1.8$	${\log }_{10}\left(\mathrm{sSFR}\right)=[-8.9,-8.4]$	15	1.54	10.13	-8.80	0.53	${1.8}_{-0.4}^{+0.4}$	${2.1}_{-0.4}^{+0.5}$	${135}_{-21}^{+13}$	${60}_{-12}^{+3}$
	$[-8.4,-7.0]$	2	1.50	10.01	-8.37	0.67	${2.3}_{-0.7}^{+4.0}$	${2.8}_{-0.9}^{+5.2}$	${199}_{-58}^{+60}$	${72}_{-30}^{+16}$
	$[6.0,8.6]$	20	1.54	9.94	-9.08	0.59	${1.6}_{-0.4}^{+0.3}$	${2.0}_{-0.5}^{+0.3}$	${108}_{-27}^{+7}$	${59}_{-10}^{+3}$
	${\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)=[8.6,9.1]$	19	1.54	10.29	-8.92	0.53	${2.3}_{-0.6}^{+0.3}$	${2.8}_{-0.8}^{+0.3}$	${165}_{-13}^{+26}$	${70}_{-11}^{+9}$
	$[9.1,13.0]$	2	1.53	10.49	-8.85	0.40	${2.3}_{-0.8}^{+3.8}$	${2.8}_{-1.2}^{+4.3}$	${162}_{-86}^{+63}$	${89}_{-21}^{+87}$
	Full redshift bin	41	1.54	10.13	-8.97	0.57	${1.9}_{-0.3}^{+0.2}$	${2.3}_{-0.4}^{+0.2}$	${135}_{-10}^{+14}$	${67}_{-11}^{+-0}$
	$[8,9.75]$	9	2.29	9.52	-8.45	0.39	${1.4}_{-0.6}^{+0.6}$	${1.8}_{-0.8}^{+0.7}$	${98}_{-38}^{+6}$	${62}_{-15}^{+9}$
	${\log }_{10}\left(M_{\ast }\right)=[9.75,10.25]$	19	2.28	10.02	-8.81	0.46	${1.0}_{-0.1}^{+0.4}$	${1.5}_{-0.3}^{+0.2}$	${94}_{-10}^{+19}$	${67}_{-8}^{+5}$
	$[10.25,12]$	22	2.27	10.51	-8.72	0.63	${1.4}_{-0.3}^{+0.1}$	${1.7}_{-0.4}^{+0.1}$	${111}_{-2}^{+30}$	${82}_{-10}^{+3}$
	$[-11,-8.9]$	10	2.21	10.21	-9.07	0.44	${1.5}_{-0.3}^{+0.5}$	${1.8}_{-0.4}^{+0.6}$	${130}_{-32}^{+7}$	${76}_{-13}^{+5}$
$2.0\le z\le 2.6$	${\log }_{10}\left(\mathrm{sSFR}\right)=[-8.9,-8.4]$	32	2.29	10.25	-8.71	0.56	${1.4}_{-0.3}^{+0.1}$	${1.7}_{-0.3}^{+0.2}$	${103}_{-5}^{+20}$	${71}_{-4}^{+7}$
	$[-8.4,-7.0]$	8	2.36	9.87	-8.18	0.48	${0.8}_{-0.2}^{+0.3}$	${1.0}_{-0.2}^{+0.4}$	${77}_{-29}^{+9}$	${62}_{-4}^{+8}$
	$[6.0,8.6]$	14	2.21	9.88	-8.78	0.55	${1.5}_{-0.4}^{+0.3}$	${1.8}_{-0.4}^{+0.5}$	${114}_{-22}^{+6}$	${67}_{-13}^{+5}$
	${\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)=[8.6,9.1]$	27	2.28	10.24	-8.72	0.54	${1.3}_{-0.3}^{+0.1}$	${1.6}_{-0.4}^{+0.1}$	${97}_{-9}^{+19}$	${67}_{-4}^{+7}$
	$[9.1,13.0]$	9	2.32	10.59	-8.59	0.41	${1.5}_{-0.5}^{+0.1}$	${1.9}_{-0.8}^{+0.1}$	${129}_{-31}^{+28}$	${88}_{-15}^{+2}$
	Full redshift bin	50	2.28	10.16	-8.71	0.52	${1.4}_{-0.3}^{+0.0}$	${1.7}_{-0.4}^{+0.0}$	${103}_{-7}^{+11}$	${70}_{-4}^{+5}$
	$[8,9.75]$	7	3.21	9.33	-7.79	0.32	${0.8}_{-0.2}^{+0.6}$	${1.1}_{-0.3}^{+0.7}$	${57}_{-8}^{+35}$	${69}_{-18}^{+3}$
	${\log }_{10}\left(M_{\ast }\right)=[9.75,10.25]$	7	3.24	9.92	-8.36	0.32	${0.8}_{-0.3}^{+0.3}$	${1.1}_{-0.3}^{+0.3}$	${85}_{-24}^{+20}$	${74}_{-6}^{+15}$
	$[10.25,12]$	0	—	—	—	—	—	—	—	—
	$[-11,-8.9]$	0	—	—	—	—	—	—	—	—
$2.9\le z\le 3.8$	${\log }_{10}\left(\mathrm{sSFR}\right)=[-8.9,-8.4]$	5	3.23	9.82	-8.47	0.32	${1.1}_{-0.3}^{+0.3}$	${1.4}_{-0.5}^{+0.4}$	${69}_{-16}^{+28}$	${69}_{-14}^{+4}$
	$[-8.4,-7.0]$	9	3.22	9.53	-7.79	0.29	${0.7}_{-0.2}^{+0.4}$	${1.0}_{-0.2}^{+0.5}$	${84}_{-29}^{+17}$	${77}_{-11}^{+13}$
	$[6.0,8.6]$	1	3.16	9.31	-8.47	0.32	${1.2}_{-0.8}^{+2.9}$	${1.6}_{-1.1}^{+3.4}$	${57}_{-47}^{+36}$	${37}_{-22}^{+12}$
	${\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)=[8.6,9.1]$	8	3.21	9.75	-8.40	0.34	${0.8}_{-0.1}^{+0.4}$	${1.0}_{-0.1}^{+0.5}$	${77}_{-15}^{+16}$	${70}_{-8}^{+5}$
	$[9.1,13.0]$	5	3.32	9.81	-7.59	0.20	${0.6}_{-0.2}^{+0.6}$	${1.1}_{-0.5}^{+0.5}$	${85}_{-32}^{+31}$	${77}_{-1}^{+29}$
	Full redshift bin	14	3.23	9.74	-8.34	0.32	${0.8}_{-0.1}^{+0.3}$	${1.1}_{-0.2}^{+0.4}$	${84}_{-20}^{+7}$	${73}_{-8}^{+6}$

Table 3. Table 3: Median kinematic properties, Unresolved/Misaligned and Undetected Rotation Galaxies

			Median Properties
Redshift	Bin	$N$	$⟨z⟩$	$⟨{\log }_{10}\left(M_{\ast }\right)⟩$	$⟨{\log }_{10}\left(\mathrm{sSFR}\right)⟩$	$⟨{\log }_{10}\left(R_E\right)⟩$	$⟨{\left(V/{\sigma }_{V,0}\right)}_{R_E}⟩$aa$V/{\sigma }_{V,0}$ for the unresolved/misaligned and velocity-limit galaxies without constrained upper uncertainties (e.g., fit slope is consistent with zero at maximum grid ${\left(V/{\sigma }_{V,0}\right)}_{R_E}=10$) have upper uncertainties marked as 99.	$⟨{\left(V/{\sigma }_{V,0}\right)}_{2.2}⟩$a,ba,bfootnotemark:
			—	$\left[{\log }_{10}\left(M_{\odot }\right)\right]$	$\left[{\log }_{10}\left({\mathrm{yr}}^{-1}\right)\right]$	$\left[{\log }_{10}\left(\mathrm{kpc}\right)\right]$	—	—
	$[8,9.75]$	58	1.53	9.50	-8.78	0.34	${2.7}_{-1.1}^{+1.8}$	${2.8}_{-1.1}^{+1.9}$
	${\log }_{10}\left(M_{\ast }\right)=[9.75,10.25]$	51	1.54	9.95	-8.95	0.40	${2.3}_{-1.1}^{+2.5}$	${2.4}_{-1.1}^{+2.6}$
	$[10.25,12]$	54	1.47	10.55	-9.13	0.58	${8.6}_{-4.6}^{+99}$	${9.0}_{-4.8}^{+99}$
	$[-11,-8.9]$	90	1.49	10.15	-9.17	0.51	${6.3}_{-2.7}^{+99}$	${6.6}_{-2.8}^{+99}$
$1.3\le z\le 1.8$	${\log }_{10}\left(\mathrm{sSFR}\right)=[-8.9,-8.4]$	60	1.53	9.68	-8.70	0.38	${2.3}_{-0.8}^{+1.2}$	${2.4}_{-0.8}^{+1.2}$
	$[-8.4,-7.0]$	13	1.57	9.83	-8.32	0.30	${1.6}_{-1.6}^{+1.1}$	${1.7}_{-1.7}^{+1.1}$
	$[6.0,8.6]$	70	1.53	9.64	-9.06	0.52	${3.1}_{-1.3}^{+4.4}$	${3.2}_{-1.4}^{+4.6}$
	${\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)=[8.6,9.1]$	71	1.50	10.23	-8.86	0.44	${3.1}_{-1.2}^{+4.2}$	${3.2}_{-1.2}^{+4.4}$
	$[9.1,13.0]$	22	1.52	10.23	-8.79	0.20	${2.2}_{-1.5}^{+1.6}$	${2.3}_{-1.6}^{+1.7}$
	Full redshift bin	163	1.52	9.93	-8.94	0.44	${3.8}_{-1.1}^{+2.2}$	${4.0}_{-1.1}^{+2.3}$
	$[8,9.75]$	99	2.24	9.53	-8.65	0.21	${1.0}_{-1.0}^{+0.9}$	${1.0}_{-1.0}^{+0.9}$
	${\log }_{10}\left(M_{\ast }\right)=[9.75,10.25]$	133	2.29	9.93	-8.76	0.31	${1.8}_{-0.8}^{+1.2}$	${1.9}_{-0.8}^{+1.2}$
	$[10.25,12]$	51	2.31	10.44	-8.76	0.47	${2.6}_{-1.7}^{+2.5}$	${2.7}_{-1.8}^{+2.6}$
	$[-11,-8.9]$	77	2.29	9.95	-9.04	0.35	${4.7}_{-2.2}^{+99}$	${4.9}_{-2.3}^{+99}$
$2.0\le z\le 2.6$	${\log }_{10}\left(\mathrm{sSFR}\right)=[-8.9,-8.4]$	147	2.28	9.85	-8.67	0.33	${1.5}_{-0.6}^{+0.7}$	${1.6}_{-0.6}^{+0.7}$
	$[-8.4,-7.0]$	59	2.28	9.81	-8.28	0.19	${0.0}_{-0.0}^{+1.1}$	${0.0}_{-0.0}^{+1.1}$
	$[6.0,8.6]$	72	2.26	9.77	-8.84	0.55	${0.4}_{-0.4}^{+1.6}$	${0.4}_{-0.4}^{+1.7}$
	${\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)=[8.6,9.1]$	122	2.29	9.85	-8.70	0.32	${1.5}_{-0.7}^{+0.7}$	${1.6}_{-0.7}^{+0.7}$
	$[9.1,13.0]$	89	2.29	9.91	-8.45	0.10	${0.3}_{-0.3}^{+0.7}$	${0.3}_{-0.3}^{+0.7}$
	Full redshift bin	283	2.29	9.86	-8.70	0.31	${1.9}_{-0.6}^{+0.6}$	${2.0}_{-0.6}^{+0.6}$
	$[8,9.75]$	55	3.23	9.53	-8.42	0.20	${1.5}_{-1.5}^{+5.2}$	${1.6}_{-1.6}^{+5.4}$
	${\log }_{10}\left(M_{\ast }\right)=[9.75,10.25]$	56	3.17	9.94	-8.65	0.29	${1.9}_{-1.0}^{+1.3}$	${2.0}_{-1.0}^{+1.4}$
	$[10.25,12]$	19	3.23	10.39	-8.68	0.35	${2.5}_{-0.4}^{+99}$	${2.6}_{-0.4}^{+99}$
	$[-11,-8.9]$	16	3.13	9.92	-8.97	0.29	${0.0}_{-0.0}^{+6.0}$	${0.0}_{-0.0}^{+6.2}$
$2.9\le z\le 3.8$	${\log }_{10}\left(\mathrm{sSFR}\right)=[-8.9,-8.4]$	71	3.19	9.91	-8.65	0.29	${2.0}_{-0.9}^{+1.4}$	${2.1}_{-0.9}^{+1.5}$
	$[-8.4,-7.0]$	43	3.30	9.66	-8.08	0.23	${0.5}_{-0.5}^{+1.4}$	${0.5}_{-0.5}^{+1.5}$
	$[6.0,8.6]$	18	3.17	9.69	-8.63	0.53	${0.0}_{-0.0}^{+0.5}$	${0.0}_{-0.0}^{+0.5}$
	${\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)=[8.6,9.1]$	58	3.19	9.90	-8.66	0.30	${1.0}_{-1.0}^{+2.6}$	${1.0}_{-1.0}^{+2.7}$
	$[9.1,13.0]$	54	3.24	9.82	-8.41	0.06	${1.4}_{-1.1}^{+1.0}$	${1.5}_{-1.1}^{+1.0}$
	Full redshift bin	130	3.21	9.84	-8.53	0.26	${1.6}_{-0.8}^{+1.0}$	${1.7}_{-0.8}^{+1.0}$

Table 4. Table 4: Median mass offsets a a Including all galaxies, where M dyn subscript 𝑀 dyn M_{\mathrm{dyn}} for unresolved/misaligned and velocity-limit galaxies are calculated using the median V / σ 𝑉 𝜎 V/\sigma for each respective bin.

			$⟨z⟩$	$⟨{\log }_{10}\left(M_{\ast }\right)⟩$	$⟨{\log }_{10}\left(\mathrm{sSFR}\right)⟩$	$⟨{\log }_{10}\left(R_E\right)⟩$	$⟨{\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)⟩$	$⟨\Delta {\log }_{10}M⟩$
Redshift	Bin	$N$	—	$\left[{\log }_{10}\left(M_{\odot }\right)\right]$	$\left[{\log }_{10}\left({\mathrm{yr}}^{-1}\right)\right]$	$\left[{\log }_{10}\left(\mathrm{kpc}\right)\right]$	$\left[{\log }_{10}\left(M_{\odot }{\mathrm{kpc}}^{-2}\right)\right]$	[dex]
	$[8,9.75]$	63	1.54	9.50	-8.81	0.33	8.40	${0.16}_{-0.03}^{+0.13}$
	${\log }_{10}\left(M_{\ast }\right)=[9.75,10.25]$	73	1.54	9.97	-8.94	0.45	8.61	${0.19}_{-0.06}^{+0.03}$
	$[10.25,12]$	68	1.48	10.56	-9.13	0.61	8.84	${0.19}_{-0.02}^{+0.07}$
	$[-11,-8.9]$	114	1.50	10.15	-9.16	0.53	8.58	${0.24}_{-0.05}^{+0.03}$
	${\log }_{10}\left(\mathrm{sSFR}\right)=[-8.9,-8.4]$	75	1.54	9.82	-8.72	0.39	8.73	${0.10}_{-0.02}^{+0.03}$
$1.3\le z\le 1.8$	$[-8.4,-7.0]$	15	1.53	9.84	-8.34	0.31	8.77	${0.15}_{-0.20}^{+0.09}$
	$R_E=[0.0,2.5]$	78	1.52	9.67	-8.81	0.29	8.82	$-{0.04}_{-0.03}^{+0.07}$
	$[2.5,15.0]$	126	1.53	10.22	-9.09	0.58	8.58	${0.26}_{-0.04}^{+0.03}$
	$[6.0,8.6]$	90	1.53	9.77	-9.06	0.54	8.34	${0.31}_{-0.05}^{+0.03}$
	${\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)=[8.6,9.1]$	90	1.52	10.23	-8.88	0.46	8.84	${0.13}_{-0.03}^{+0.03}$
	$[9.1,13.0]$	24	1.52	10.31	-8.79	0.20	9.22	$-{0.02}_{-0.08}^{+0.09}$
	$[8,9.75]$	108	2.24	9.53	-8.62	0.23	8.72	$-{0.07}_{-0.02}^{+0.03}$
	${\log }_{10}\left(M_{\ast }\right)=[9.75,10.25]$	152	2.29	9.94	-8.76	0.34	8.90	$-{0.11}_{-0.02}^{+0.02}$
	$[10.25,12]$	73	2.30	10.46	-8.76	0.54	9.00	$-{0.03}_{-0.05}^{+0.05}$
	$[-11,-8.9]$	87	2.29	10.00	-9.04	0.38	8.74	${0.11}_{-0.10}^{+0.04}$
	${\log }_{10}\left(\mathrm{sSFR}\right)=[-8.9,-8.4]$	179	2.28	9.89	-8.68	0.35	8.80	$-{0.08}_{-0.03}^{+0.03}$
$2.0\le z\le 2.6$	$[-8.4,-7.0]$	67	2.29	9.81	-8.26	0.23	9.18	$-{0.21}_{-0.04}^{+0.02}$
	$R_E=[0.0,2.5]$	198	2.27	9.78	-8.65	0.22	9.06	$-{0.12}_{-0.02}^{+0.03}$
	$[2.5,15.0]$	135	2.29	10.12	-8.79	0.57	8.63	$-{0.03}_{-0.06}^{+0.03}$
	$[6.0,8.6]$	86	2.25	9.80	-8.83	0.55	8.39	${0.06}_{-0.06}^{+0.06}$
	${\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)=[8.6,9.1]$	149	2.29	9.92	-8.71	0.35	8.83	$-{0.09}_{-0.02}^{+0.04}$
	$[9.1,13.0]$	98	2.29	9.94	-8.46	0.14	9.31	$-{0.26}_{-0.03}^{+0.02}$
	$[8,9.75]$	62	3.23	9.52	-8.39	0.22	9.01	$-{0.21}_{-0.08}^{+0.09}$
	${\log }_{10}\left(M_{\ast }\right)=[9.75,10.25]$	63	3.18	9.94	-8.58	0.29	9.02	$-{0.10}_{-0.06}^{+0.02}$
	$[10.25,12]$	19	3.23	10.39	-8.68	0.35	9.35	$-{0.25}_{-0.13}^{+0.17}$
	$[-11,-8.9]$	16	3.13	9.92	-8.97	0.29	8.93	${0.01}_{-0.07}^{+0.06}$
	${\log }_{10}\left(\mathrm{sSFR}\right)=[-8.9,-8.4]$	76	3.20	9.91	-8.64	0.29	8.91	$-{0.11}_{-0.05}^{+0.02}$
$2.9\le z\le 3.8$	$[-8.4,-7.0]$	52	3.27	9.65	-8.08	0.24	9.28	$-{0.28}_{-0.10}^{+0.03}$
	$R_E=[0.0,2.5]$	107	3.22	9.76	-8.46	0.18	9.18	$-{0.23}_{-0.04}^{+0.04}$
	$[2.5,15.0]$	37	3.21	9.96	-8.62	0.49	8.69	$-{0.05}_{-0.03}^{+0.12}$
	$[6.0,8.6]$	19	3.16	9.66	-8.62	0.53	8.46	$-{0.04}_{-0.05}^{+0.01}$
	${\log }_{10}\left({\Sigma }_{\mathrm{bar}}\right)=[8.6,9.1]$	66	3.20	9.88	-8.59	0.31	8.87	$-{0.07}_{-0.10}^{+0.04}$
	$[9.1,13.0]$	59	3.24	9.81	-8.38	0.06	9.42	$-{0.29}_{-0.08}^{+0.03}$

Equations

$$V_{\mathrm{circ}}(r)=\sqrt{V(r{)}^2+2\left(r/r_s\right){\sigma }_{V,0}^2},$$

$$M_{\mathrm{dyn}}=k_{\mathrm{tot}}(R_E)\frac{V_{\mathrm{circ}}(R_E{)}^2{R}_E}{G},$$

$$V_{\mathrm{circ}}(R_E){|}_{(V/{\sigma }_{V,0})}=\sqrt{\left[\frac{(V/{\sigma }_{V,0}{)}_{R_E}^2+3.35}{(V/{\sigma }_{V,0}{)}_{R_E}^2+1}\right]}{V}_{\mathrm{rms}}(R_E{)}_{1D,\mathrm{corr}},$$

Full text

Galaxy kinematics (602); Galaxy dynamics (591); High-redshift galaxies (734); Galaxy evolution (594)

1 Introduction

A key open question in galaxy formation and evolution is how galaxy structures arise and evolve over time. Today’s massive star-forming galaxies are assumed to form by the collapse of baryons within dark matter halos (e.g., White & Rees 1978, Fall & Efstathiou 1980, Blumenthal et al. 1984), resulting in thin, smooth stellar disks. However, the exact details of how baryons and dark matter interact throughout the galaxy disk formation process are not well understood. Constraining these physical processes and testing different formation models require direct studies of galaxies at earlier times.

Recent work shows that massive star-forming galaxies in the early universe look very different from their local counterparts. In particular, at $z\sim 1-3$, the peak of cosmic star formation rate (SFR) density in the universe (e.g., Madau & Dickinson 2014), massive star-forming galaxies are generally smaller (e.g., Williams et al. 2010, van der Wel et al. 2014a), have large clumps (e.g., Elmegreen et al. 2007, Law et al. 2007, Genzel et al. 2011, Förster Schreiber et al. 2011, Guo et al. 2015, 2018), have high gas fractions (e.g., Daddi et al. 2008, Tacconi et al. 2008, 2018), and are baryon dominated within their disk extent (e.g., van Dokkum et al. 2015, Wuyts et al. 2016, Genzel et al. 2017, Lang et al. 2017, Übler et al. 2017, Übler et al. 2018). Despite these differences, high-redshift massive galaxies do appear to have rotating gas disks, though they are thicker (e.g., Elmegreen & Elmegreen 2006) and have higher intrinsic velocity dispersions and increased turbulence compared to local massive star-forming galaxies (e.g., Förster Schreiber et al. 2006, 2009, Epinat et al. 2008, Newman et al. 2013, Green et al. 2014, Wisnioski et al. 2015, Simons et al. 2016).

Current theoretical models suggest that the thick, gas-rich disks of $z\sim 1-3$ massive star-forming galaxies are assembled through smooth, cold-mode accretion or minor mergers (e.g., Kereš et al. 2005, 2009, Dekel & Birnboim 2006, Davé 2008, Dekel et al. 2009, Oser et al. 2010, Cacciato et al. 2012, Ceverino et al. 2012, Danovich et al. 2015, Rodriguez-Gomez et al. 2015, Correa et al. 2018). In this framework, the high turbulence and clumpy morphologies of these galaxies at $z\sim 1-3$ could reflect higher average gas fractions compared to local galaxies (e.g., Dekel et al. 2009, Bournaud et al. 2011, Genzel et al. 2011, Genel et al. 2012).

Recent instrumentation advances, including multiplexing near-infrared (NIR) spectrographs such as MOSFIRE (McLean et al. 2010, 2012) and KMOS (Sharples et al. 2004, 2013), have enabled large kinematic studies at high redshifts to test these theoretical models. Surveys with MOSFIRE and KMOS (e.g., MOSDEF, Kriek et al. 2015; SIGMA, Simons et al. 2016; ZFIRE, Nanayakkara et al. 2016; KMOS${}^{\hbox{\scriptsize{3D}}}$, Wisnioski et al. 2015; KROSS, Stott et al. 2016; and KDS, Turner et al. 2017) now provide the kinematics for thousands of galaxies at $z\sim 1-3$, augmenting more detailed studies of smaller samples with higher, adaptive-optics-assisted spatial resolution (e.g., SINS/zC-SINF, Förster Schreiber et al. 2009, 2018; MASSIV, e.g., Contini et al. 2012; AMAZE/LSD, e.g., Gnerucci et al. 2011). While these surveys have greatly expanded our understanding of high-redshift galaxies, many challenges to interpreting these results remain.

First, the majority of previous high-redshift kinematic studies are conducted using ground-based, seeing-limited instruments (e.g., MOSFIRE and KMOS). The low spatial resolution of these observations can mask rotation signatures in small and lower-mass galaxies (Newman et al., 2013), impacting most galaxies with $M_{*}\lesssim 10^{10}M_{\odot}$ at $z\sim 1-3$ as these galaxies are unresolved under seeing-limited conditions (van der Wel et al. 2014a). Second, multiobject slit spectrographs (e.g., MOSFIRE) often have a constant position angle for all slits in a mask, resulting in random galaxy orientations within the slits that can further mask rotation signals. Third, fully constraining star-forming galaxy formation models requires observations of lower-mass galaxies at high redshifts. The progenitors of today’s massive star-forming disk galaxies have masses $M_{*}\sim 10^{9}-10^{10}M_{\odot}$ at $z\sim 2$ (e.g., Leja et al. 2013, van Dokkum et al. 2013), but currently kinematic observations of these early, low-mass galaxies are limited to very small sample sizes.

Fortunately, high-resolution space-based imaging allows us to constrain the structures of high-redshift galaxies by combining the information from the detailed imaging with seeing-limited spectra (e.g., Price et al. 2016). By leveraging morphology from space-based imaging, it is possible to account for how much of the galaxy falls outside of the slit due to misalignment and for the effects of seeing blurring, which enables us to better understand the galaxy components and kinematics that are captured in the seeing-limited spectra. This approach can be applied to galaxies both with and without spatially resolved kinematics and can thus be used to study the internal dynamics for galaxies over a wide range of masses and sizes at high redshifts, as in Price et al. 2016.

In this paper, we use observations from the MOSFIRE Deep Evolution Field (MOSDEF) survey (Kriek et al., 2015) to study the dynamical and baryonic masses and internal kinematic structures of a sample of 681 star-forming galaxies at $z\sim 1.4-3.8$. Rotation is robustly detected in 105 galaxies, after we initially restrict our analysis to galaxies that are spatially resolved along the slit and that have rough alignment between the galaxy major axis and the slit. We use detailed structural information from Hubble Space Telescope (HST) imaging from the CANDELS survey (Grogin et al. 2011, Koekemoer et al. 2011) and other ancillary information to statistically constrain the kinematics of the 95 spatially resolved, slit-aligned galaxies without detected rotation and of the remaining 481 spatially unresolved or slit-misaligned galaxies. We investigate how the derived intrinsic rotation velocities ($V$), intrinsic velocity dispersions ($\sigma_{V,0}$), and ratio of intrinsic ordered to unordered motion ($V/\sigma_{V,0}$) correlate with other galaxy properties, and how these kinematic properties change over cosmic time. We also compare the derived dynamical masses with the galaxies’ baryonic masses to infer the evolution of the baryon and dark matter fraction in galaxies as a function of galaxy properties and redshift.

We adopt a $\Lambda$CDM cosmology with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and $H_{0}=70\mathrm{\,km}\mathrm{\,s^{-1}}\mathrm{\,Mpc^{-1}}$ throughout this work.

2 Data

2.1 The MOSDEF Survey

This work is based on the complete data set from the MOSDEF survey (Kriek et al., 2015), which was carried out from 2012 December to 2016 May using the MOSFIRE spectrograph (McLean et al., 2012) on the 10 m Keck I telescope. In total, the survey obtained rest-frame optical, moderate-resolution ($R=3000-3600$) spectra for $\,\sim\!1500$ $H$-band-selected galaxies in the CANDELS (Grogin et al. 2011, Koekemoer et al. 2011) fields. The full survey details, including targeting, observational strategy, data reduction, success rate, sensitivities, redshift measurements, and other sample properties, are given in Kriek et al. (2015).

Structural parameters for the MOSDEF galaxies are measured by van der Wel et al. (2014a) from HST/F160W images from the CANDELS survey using Galfit (Peng et al., 2010). These parameters include the effective radius, $R_{E}$ (using the semi-major axis length as $R_{E}$), the Sérsic index, $n$ (Sérsic, 1968), the axis ratio, $b/a$, and the major-axis position angle.

The stellar masses and other stellar population parameters, including spectral energy distribution (SED) SFRs, are determined by fitting the $0.3-8.0\mu\mathrm{m}$ multiwavelength photometry from the 3D-HST survey (Brammer et al. 2012, Skelton et al. 2014, Momcheva et al. 2016) with stellar population models, in the same way as described in Price et al. (2016) (see also Kriek et al. 2015). Briefly, FAST (Kriek et al., 2009) is used to fit the flexible stellar population models (Conroy et al. 2009, Conroy & Gunn 2010) to the multiwavelength photometry, using the MOSFIRE redshifts and assuming a Chabrier (2003) stellar initial mass function (IMF), a Calzetti et al. (2000) dust attenuation curve, delayed exponentially declining star formation histories, and solar metallicity. We additionally correct the stellar masses for differences between the Galfit and total photometric magnitudes following Taylor et al. (2010) (as in Equation (1) of Price et al. 2016), ensuring self-consistency between the radii and stellar mass measurements.

Emission line fluxes (e.g., $\textrm{H}\alpha$, [Oiii], $\textrm{H}\beta$) are calculated from the optimally extracted MOSFIRE 1D spectra by fitting the lines and underlying continuum simultaneously with Gaussian profiles and a linear component. The $\textrm{H}\alpha$ and $\textrm{H}\beta$ fluxes are additionally corrected for the underlying Balmer absorption, using the best-fit stellar population models. More details on emission line measurements are given in Kriek et al. (2015) and Reddy et al. (2015).

We use the following ranked approach to measure SFRs and estimate gas masses ($M_{\mathrm{gas}}$) for our sample. First, if both $\textrm{H}\alpha$ and $\textrm{H}\beta$ are detected ($\mathrm{S/N}\geq 3$), the spectrum transmission at the lines is at least 50% of the maximum transmission, and both lines have reliable fluxes (i.e., removing lines at the edge of the spectrum with unreliable fluxes), then Balmer-absorption-corrected $\textrm{H}\alpha$ fluxes are corrected for dust attenuation using the Balmer decrement (320 galaxies, median $f_{\textrm{H}\alpha}/f_{\textrm{H}\beta}\sim 3.43$), assuming a Cardelli et al. (1989) extinction curve (see Reddy et al. 2015). $\textrm{H}\alpha$ luminosities are then calculated from the dust-corrected $\textrm{H}\alpha$ fluxes, and $\textrm{H}\alpha$ SFRs are calculated from the luminosities using the relation of Hao et al. (2011) for a Chabrier (2003) IMF (Shivaei et al., 2015).

Next, if only one Balmer line ($\textrm{H}\alpha$ or $\textrm{H}\beta$) is detected with signal-to-noise ratio $\mathrm{S/N}\geq 3$ and has a reliable flux (283 galaxies), then we assume that the dust attenuation of the nebular regions is related to the continuum attenuation using the relation between $A_{V,\mathrm{neb,Calzetti}}$ and $A_{V,\mathrm{cont,Calzetti}}$ by Price et al. (2014). As we have chosen to adopt the Cardelli et al. (1989) curve for the nebular attenuation, we convert this relation to account for a Cardelli et al. (1989) nebular and a Calzetti et al. (2000) continuum attenuation curve.111As we determine the color excess from the observed Balmer decrement using $E(\mathrm{B-V})_{\mathrm{neb}}=\log_{10}\left[(\textrm{H}\alpha/\textrm{H}\beta)_{\mathrm{obs}}/2.86\right]/\left[0.4\,(k(\lambda_{\textrm{H}\beta})-k(\lambda_{\textrm{H}\alpha}))\right]$ for some reddening curve $k(\lambda)$, for the same observed Balmer decrement we have $E(\mathrm{B-V})_{\mathrm{neb,Cardelli}}=1.18E(\mathrm{B-V})_{\mathrm{neb,Calzetti}}$. Price et al. (2014) Eq. 2 is equivalent to $E(\mathrm{B-V})_{\mathrm{neb}}=1.86\,E(\mathrm{B-V})_{\mathrm{cont}}$ as the Calzetti et al. (2000) curve was adopted for both the nebular emission and stellar continuum, so this converts to $E(\mathrm{B-V})_{\mathrm{neb,Cardelli}}=2.2\,E(\mathrm{B-V})_{\mathrm{cont,Calzetti}}$. Since $A_{\mathrm{V}}=R_{\mathrm{V}}\times E(\mathrm{B-V})$, as $R_{\mathrm{V}}=3.1$ for the Cardelli et al. (1989) curve and $R_{\mathrm{V}}=4.05$ for the Calzetti et al. (2000) curve, the final converted relation is $A_{V,\mathrm{neb,Cardelli}}=1.68\,A_{V,\mathrm{cont,Calzetti}}$. We then use the resulting attenuation $A_{V,\mathrm{neb,Cardelli}}$ to correct the absorption-corrected Balmer line flux. The Balmer luminosities are then calculated from the dust-corrected line fluxes. If only $\textrm{H}\beta$ is detected, we convert the dust-corrected $\textrm{H}\beta$ luminosity to an $\textrm{H}\alpha$ luminosity assuming $L_{\textrm{H}\alpha}/L_{\textrm{H}\beta}=2.86$ (Osterbrock & Ferland, 2006). As with the above, $\textrm{H}\alpha$ (or $\textrm{H}\beta$) SFRs are then determined from the dust-corrected Balmer luminosities using the Hao et al. (2011) relation.

Finally, if neither Balmer line is detected (78 galaxies), SFRs from SED fitting are adopted. In the absence of detected Balmer lines, these SED SFRs should be fairly reasonable to use, as Shivaei et al. (2016) show that there is general consistency between Balmer decrement-corrected $\textrm{H}\alpha$ and SED SFRs for a subset of the MOSDEF sample. However, we do note that Reddy et al. (2015) find $\textrm{H}\alpha$ SFRs are systematically higher than SED SFRs for galaxies with high SFRs or specific SFRs (sSFRs) for a separate MOSDEF sample subset. To understand the impact of assuming different SFR indicators in this ‘ladder’ method, we examine the agreement between the SFR measurements for the galaxies at $z\sim 2$ that have detections of both $\textrm{H}\alpha$ and $\textrm{H}\beta$ in Appendix B. We find that the indicators are generally in good agreement, though the SED SFRs are $\sim-0.23\mathrm{\,dex}$ lower than the Balmer-decrement $\textrm{H}\alpha$ SFRs. (We discuss further implications of these SFR offsets in Section 5.2 and Appendix B.) While the SFR indicators are not perfectly matched, this staggered method nonetheless allows us to expand the sample selection from Price et al. (2016) to include galaxies over a wider range of redshifts, as $\textrm{H}\alpha$ measurements at $z\sim 3$ are not accessible from the ground.

We then estimate gas masses for every galaxy using the relation $\Sigma_{\mathrm{SFR}}\propto\Sigma_{\mathrm{gas}}^{N}$, with $N=1.4$ following Kennicutt (1998), as this slope is intermediate between previously measured Schmidt-law slopes at $z\sim 1-3$ (see discussion in Section 5.1.1). Alternatively, gas masses could be estimated using the scaling relations by Tacconi et al. (2018). However, we choose to use the Kennicutt relation, as we primarily adopt $\textrm{H}\alpha$ SFRs here and not UV+IR SFRs (though we discuss in Section 5.1.2 how adopting gas masses from these scaling relations would impact our analysis). Here we use $\Sigma_{\mathrm{SFR}}=\mathrm{SFR}/(2\pi R_{E}^{2})$ and $\Sigma_{\mathrm{gas}}=M_{\mathrm{gas}}/(2\pi R_{E}^{2})$,222We erroneously used $\Sigma_{\mathrm{SFR}}=\mathrm{SFR}/(\pi R_{E}^{2})$ and $\Sigma_{\mathrm{gas}}=M_{\mathrm{gas}}/(\pi R_{E}^{2})$ in Price et al. (2016). Using the correct surface density definitions and the Kennicutt (1998) relation for a Chabrier IMF, the median $\log_{10}(M_{\mathrm{dyn}}/M_{\mathrm{baryon}})$ for our previous work should be $\Delta\log_{10}M=-0.02$ (lower by $-0.06\mathrm{\,dex}$), implying no dark matter (instead of $f_{\mathrm{DM}}=8\%$). Similarly, the intercept of the $M_{\mathrm{baryon}}$-$S_{0.5}$ relation should be lower by 0.01 dex in $\log_{10}(S_{0.5})$. Nonetheless, our primary result that galaxies at $z\sim 1.5-2$ are very baryon dominated within their effective radii remains unchanged, as does the conclusion that a Chabrier IMF is more consistent with the dynamical masses than a Salpeter IMF.  where we adopt the best-available SFR and $R_{E}$ as the best-fit Galfit semi-major axis, assuming that the emission-line region coincides with the stellar continuum.

2.2 Sample Selection

We select a kinematics sample from the full MOSDEF survey using the following criteria. Galaxies are selected in the redshift ranges $1.34\leq z\leq 1.75$, $2.075\leq z\leq 2.6$, and $2.9\leq z\leq 3.8$, to include galaxies in all three redshift ranges covered by the survey, with the best possible $\textrm{H}\alpha$ coverage for the lower two ranges. We next require that each galaxy has at least one emission line of $\textrm{H}\alpha$, $\textsc{[Oiii]}5007\Angstrom$, or $\textrm{H}\beta$ detected with $\mathrm{S/N}\geq 3$ and without severe skyline contamination. The sample is further restricted to only galaxies with HST/F160W coverage, as we require structural parameter measurements. Where possible, we use $\textrm{H}\alpha$ to measure kinematics to provide a fair comparison with previous studies using $\textrm{H}\alpha$ kinematics (e.g., Förster Schreiber et al. 2009, Wisnioski et al. 2015, Stott et al. 2016), but if $\textrm{H}\alpha$ fails the S/N or skyline contamination cuts, we include galaxies with clean $\textsc{[Oiii]}5007\Angstrom$ or $\textrm{H}\beta$ detections (in order of preference). We prioritize $\textsc{[Oiii]}5007\Angstrom$ over $\textrm{H}\beta$ to provide the most direct comparison with previous $z\sim 3$ studies that use [Oiii] as the kinematic tracer (e.g., Gnerucci et al. 2011, Turner et al. 2017).

Additional selection criteria are applied to restrict the sample to galaxies with high-quality spectra and structural and stellar population parameters. First, we only include primary MOSDEF targets and exclude any serendipitously detected galaxies. Second, we only include galaxies with secure redshifts ($\mathtt{ZQUALFLAG}\geq 3$ in the final MOSDEF redshift catalog). Third, we impose stellar population and structural parameter quality cuts to ensure that the models reasonably match the data. For the SED fitting, we exclude objects with best-fit reduced chi-square $\chi^{2}_{\mathrm{red}}>10$. We then exclude galaxies with structural parameter fits where (a) Galfit failed or reached a parameter limit or (b) the difference between the Galfit and HST/F160W total magnitudes is greater than 0.5 mag. Fourth, we exclude active galactic nuclei (AGNs) from our sample based on X-ray luminosity, IRAC color, and rest-frame optical emission line ratios (Coil et al. 2015, Azadi et al. 2017), as kinematics of AGNs may include contributions from nonvirial motions in the nuclear region. Fifth, we exclude objects within the quiescent region in the UVJ diagram (Wuyts et al. 2007, Williams et al. 2009), to avoid objects where line emission comes from otherwise unidentified nuclear activity. Sixth, we remove any galaxy that appears to have an interacting counterpart at a similar redshift (i.e., objects with secondary objects within 1 arcsec, $|\Delta z|<0.1$, and with $\Delta m_{\mathrm{F160W}}\leq 1$, as well as objects with disturbed/suspicious line or spatial profiles, based on visual inspection), as these systems may have disturbed kinematics (see, e.g., Shapiro et al. 2008). Finally, any object that has a line width smaller than the instrumental resolution is excluded from this analysis.

The final kinematic sample includes 681 galaxies. Of these galaxies, kinematics are measured using $\textrm{H}\alpha$, $\textsc{[Oiii]}5007\Angstrom$, and $\textrm{H}\beta$ for 481, 195, and 5 galaxies, respectively (with some objects having two or more observations of the best-available line from different masks). Within our sample, both $\textrm{H}\alpha$ and $\textrm{H}\beta$ are detected and have reliable fluxes in 320 galaxies (i.e., Balmer-decrement -orrected $\textrm{H}\alpha$ SFRs are adopted), one Balmer line is detected with reliable flux for 283 galaxies (i.e., $\textrm{H}\alpha$ or $\textrm{H}\beta$ SFRs are estimated using $A_{V}$ from SED fitting), and no Balmer line is detected for 78 galaxies (i.e., SED SFRs are adopted).

The stellar masses, SFRs, and effective radii of our final kinematics sample are shown in Figure 1, colored by redshift range. For comparison, we show the best-fit stellar mass-SFR relation at $z\sim 1.4-2.6$ by Shivaei et al. (2015) (left panel) and the best-fit size-stellar mass relations for star-forming galaxies at $z=1.75$, $z=2.25$, and $z=2.75$ by van der Wel et al. (2014a) (right panel). The low- and medium-redshift galaxies in our sample are in excellent agreement with the best-fit relation of Shivaei et al. (2015), and overall the sample shows the expected trend of higher SFRs at higher redshifts. The galaxies also generally follow the size-mass relations measured at similar redshifts, though the galaxies at $z\sim 3$ are generally smaller than the relation measured at $z=2.75$, as our highest sample redshift range ($2.9\leq{}z\leq 3.8$) probes higher redshifts than presented by van der Wel et al. (2014a).

2.3 Sample Spatial Resolution

We next measure the point-spread function (PSF) corrected spatial extent of the highest-S/N emission line for each galaxy, as 2D kinematic fitting will only yield meaningful constraints if the emission line is resolved. By measuring the intrinsic emission line sizes from the spectra, we directly account both for size variations between a galaxy’s continuum and emission and for projection and seeing effects that impact how much of a galaxy’s light falls within the slit. Intrinsic emission line half-light diameters are measured following a similar procedure to that used in Simons et al. (2016).

First, we mask skylines and low-S/N cross-dispersion spatial rows within the 2D spectrum of each galaxy. Any continuum is then subtracted from the 2D emission line. We measure the continuum slope from a weighted linear fit to the 1D spectrum and then fit the continuum row by row using weighted linear fits where only the intercept is variable. The observed emission line FWHM is measured from a weighted fit to a Gaussian profile. The intrinsic half-light diameter, $D_{\mathrm{50,emis}}$, is then determined by subtracting the PSF FWHM in quadrature from the observed line FWHM. For galaxies with observed FWHMs or unmasked spatial ranges less than the PSF size, we set $D_{\mathrm{50,emis}}=0$.333Note that $D_{\mathrm{50,emis}}$ is only used for the spatial resolution classification and not for analysis, so for simplicity we do not adopt an upper limit. The distribution of $D_{\mathrm{50,emis}}/\mathrm{FWHM_{PSF}}$ versus the galaxy photometric major axis-slit misalignment ($|\Delta\mathrm{PA}|$) for our sample is shown in Figure 2a.

Galaxies are then classified as spatially resolved if they satisfy $D_{\mathrm{50,emis}}\geq 0.8\times\mathrm{FWHM_{PSF}}$. Simons et al. (2016) demonstrate that kinematics can be accurately recovered down to this limit for the typical central $\mathrm{S/N}\sim\!15$ per pixel for their sample, and the MOSDEF observations have similar S/N values. Their finding is consistent with the results of our kinematic recovery tests presented in Appendix A (see Figure 10), reinforcing the adoption of this resolution criterion.

3 Kinematic Measurements

The kinematic properties and dynamical masses of our sample are measured from emission lines ($\textrm{H}\alpha$, [Oiii], or $\textrm{H}\beta$), together with the structural parameters derived from the HST/F160W imaging, using the 3D models and methods developed in Price et al. (2016). In this analysis, we measure and analyze both “resolved” and “unresolved” kinematics, as discussed below. For clarity, definitions of key kinematic variables are listed in Table 1.

3.1 Resolved Kinematics

For spatially resolved galaxies, the rotation velocity and velocity dispersion can be measured from 2D spectra. However, the MOSDEF galaxies were observed with random misalignments between the galaxy major axis and the slit, resulting in reduced rotation signatures since part of the galaxy falls outside of the slit. We thus use the 3D models and methods developed by Price et al. (2016) to fit the rotation and velocity dispersion from the 2D MOSFIRE spectra, incorporating the spatial information from the HST/F160W imaging to account for the degree of slit-axis misalignment.

Full details of the kinematic models and the fitting code (MISFIT; MISaligned-slit kinematic FITting) are given in Appendix A of Price et al. (2016). In brief, the models have three free kinematic parameters: the asymptotic velocity ($V_{a}$) and turnover radius ($r_{t}$) of an arctan rotation curve model, and a constant intrinsic velocity dispersion ($\sigma_{V,0}$). The models also account for inclination (assuming an intrinsic disk axis ratio $(b/a)_{0}=0.19$, following Miller et al. 2011) and the galaxy sizes, brightness profiles, slit misalignments and kinematic aperture losses, seeing conditions, and instrumental resolution. Specifically, each model is calculated by first constructing a galaxy profile from the intrinsic kinematic and structural parameters, convolving with the PSF and instrument line-spread function, and finally extracting the kinematics within a misaligned slit. As the MISFIT models are constructed by forward modeling, the impact of beam smearing on the recovered kinematics is explicitly included. Thus, the fitting is able to directly constrain uncertainties from the range of intrinsic kinematics that match the observed data after applying beam smearing and other observational effects.

We subtract the continuum emission and mask the low-S/N rows and columns of the 2D emission lines similar to the procedure of Price et al. (2016), while adopting appropriate wavelength ranges that exclude neighboring features for the different emission lines (i.e., $\textrm{H}\alpha$, [Oiii], $\textrm{H}\beta$). Any emission lines from serendipitous objects are masked, if they fall within the fitting regions. The masked, continuum-subtracted 2D emission lines are then fit to the models using the python Markov Chain Monte Carlo (MCMC) package emcee (Foreman-Mackey et al., 2013). The MCMC sampling was conducted for 100 steps with 2000 walkers, with a burn-in phase of 50 steps. For every object, the final acceptance fraction is between 0.2 and 0.5 and the chain was run for longer than 10 times the maximum estimated parameter autocorrelation time. We then marginalize over the posterior distributions to determine the best-fit values and confidence intervals for the intrinsic $V(R_{E})$, $V_{2.2}=V(r=2.2r_{s})$ (where $r_{s}=R_{E}/1.676$ is the scale length of an exponential disk), and $\sigma_{V,0}$.

While the MISFIT 3D models directly account for the misalignment between the galaxy kinematic major axis and the slit, $\Delta\mathrm{PA}$, we lack independent measurements of the kinematic major-axis position angles. We therefore assume that the photometric and kinematic major axes are aligned. This assumption is generally reasonable, as Wisnioski et al. (2015) find that these axes are within $\Delta\mathrm{PA}\leq 30^{\circ}$ for most of their sample. Moreover, even with direct modeling of slit misalignments, objects that are very misaligned will suffer from degeneracies between the blended rotation signatures and the intrinsic velocity dispersion. In Appendix A, we use a suite of simulated galaxies and find that kinematics are recovered less accurately and have larger scatter if $|\Delta\mathrm{PA}|>45^{\circ}$ (Figure 10, fourth column). Thus, we restrict the sample for which we fit the 2D kinematics to only galaxies that are spatially resolved ($D_{\mathrm{50,emis}}\geq 0.8\times\mathrm{FWHM_{PSF}}$; Section 2.3) and relatively aligned with the slit ($|\Delta\mathrm{PA}|\leq 45^{\circ}$). The spatial resolution and slit misalignment classification of our sample is shown in Figure 2a.

Using the 2D fitting results, we then determine which spatially resolved/aligned galaxies have robust rotation detections. A galaxy is classified as having detected rotation if the best-fit rotation velocity at the effective radius, $V(R_{E})$, is (a) nonzero with 68% confidence ($|V(R_{E})|/\langle err_{V(R_{E})}\rangle\geq 1$) and    (b)    larger     than    the     instrument    resolution     limit ($2V(R_{E})/\sigma_{\mathrm{inst.\ res}}\geq 1$). We show these rotation resolution criteria for the spatially resolved/aligned galaxies in Figure 2b.

Based on these sample classification criteria, we find that 200 galaxies are spatially resolved and aligned with the slit (“resolved/aligned”). Of the resolved/aligned galaxies, rotation is detected for 105 (“resolved/aligned, rotation”), and we place limits on the rotation for the other 95 galaxies (“resolved/aligned, rotation-limit”). The remaining 481 galaxies are classified as “unresolved/misaligned”. For the galaxies with duplicate observations of the best-available line, we perform kinematic fits (either resolved or unresolved; Section 3.2) on the separate observations and then classify the object following its best-case kinematic observation(s) (resolved rotation, then rotation-limit, then unresolved/misaligned). The kinematics used in the rest of the analysis is taken to be the average of any duplicate best-case kinematics observations, weighted by the emission line S/N, or the single best-case observation, if the other observation(s) are less spatially/kinematically resolved.

The distribution of the galaxies within these three categories over redshift, stellar mass, sSFR, effective radius ($R_{E}$), and baryon mass surface density ($\Sigma_{\mathrm{bar},e}$) is shown in Figure 3. Overall, resolved/aligned galaxies with detected rotation tend to have slightly higher masses and effective radii than those galaxies with rotation limits or that are unresolved/misaligned, but all categories have similar sSFRs. Additionally, we detect rotation in relatively few galaxies at $z\sim 3$ compared to the galaxies at $z\sim 1.5$ and $z\sim 2$. This dearth of rotation-detected objects at $z\sim 3$ is likely driven by evolution of galaxy half-light sizes, resulting in fewer resolution elements at higher redshifts on average (e.g., see Fig 3c).

3.2 Unresolved Kinematics

Kinematics for galaxies classified as “unresolved/misaligned” and “resolved/aligned, rotation-limit” (Section 3.1) are measured from the integrated 1D spectra. However, while these two categories of galaxies have no detected rotation, we might expect that they have at least partial rotational support, which has likely been masked by the seeing-limited conditions (e.g., Newman et al. 2013), together with the misalignment of the galaxies and slits. Therefore, as in Price et al. (2016), we model the integrated, unresolved kinematics using the same 3D models while assuming a fixed value of $(V/\sigma_{V,0})_{R_{E}}=V(R_{E})/\sigma_{V,0}$.

Briefly, we first measure the integrated velocity dispersions $\sigma_{V,\mathrm{1D,obs}}$ by fitting the emission line ($\textrm{H}\alpha$, [Oiii], or $\textrm{H}\beta$), any neighboring lines, and the continuum simultaneously using Gaussian profiles and a linear component, assuming that the emission lines have the same width in velocity space. The measured integrated velocity dispersion $(\sigma_{V,obs})$ is then corrected for instrumental resolution by subtracting $\sigma_{V,\mathrm{sky}}=\sigma_{{inst.\ res}}$ in quadrature from $\sigma_{V,\mathrm{1D,obs}}$, which is measured from the skylines in each spectrum. Errors on $\sigma_{V,\,\mathrm{obs}}$ are estimated by repeating the fitting and correction procedure on 500 random, error-perturbed copies of each spectrum.

We then use the same MISFIT 3D kinematic models to convert the observed velocity dispersions $\sigma_{V,\mathrm{obs}}$ to an intrinsic rms velocity $V_{\mathrm{rms}}(R_{E})=\sqrt{V(R_{E})^{2}+\sigma_{V,0}^{2}}$, as described in Appendix B of Price et al. (2016). In summary, each galaxy is modeled as an inclined disk using the best-fit Galfit structural parameters and is offset from the slit by the measured $\Delta\mathrm{PA}$. Rotation and velocity dispersion are included by assuming a fixed ratio of $(V/\sigma_{V,0})_{R_{E}}$ (see Section 3.3.2) and assuming an arctan velocity curve where $r_{t}=0.4\,r_{s}=0.4\,(R_{E}/1.676)$ (Miller et al., 2011). We then compute the luminosity-weighted, seeing-convolved integrated velocity dispersion ($\sigma_{V,\mathrm{model}}$) and rms velocity at $R_{E}$ ($V_{\mathrm{rms}}(R_{E})_{\mathrm{model}}$) from the galaxy model. Finally, we calculate the composite $V_{\mathrm{rms}}(R_{E})_{\mathrm{1D,corr}}$ as in Equation (2) of Price et al. (2016), using the measured instrument resolution corrected integrated velocity dispersion ($\sigma_{V,\mathrm{obs}}$) together with the model-integrated velocity dispersion ($\sigma_{V,\mathrm{model}}$) and the model rms velocity ($V_{\mathrm{rms}}(R_{E})_{\mathrm{model}}$).

3.3 Measuring Dynamical Masses and $V/\sigma_{V,0}$

We determine the dynamical masses ($M_{\mathrm{dyn}}$) for the galaxies in our sample by combining their measured kinematics and structural information. Here we give the calculation of $M_{\mathrm{dyn}}$ from both 2D and 1D spectra, detail how $V/\sigma_{V,0}$ is constrained for galaxies without resolved and detected rotation, and compare the dynamical masses measured using both methods.

3.3.1 Resolved Rotation from 2D Spectra

For galaxies with resolved and detected rotation measured from 2D spectra, the dynamical masses are determined from their inferred rotation velocities. Furthermore, we apply an asymmetric drift correction to account for the nonnegligible pressure support in these galaxies (e.g., Epinat et al. 2009, Newman et al. 2013, Wuyts et al. 2016). The pressure-corrected circular velocity is

[TABLE]

where $r_{s}$ is the disk scale length (e.g., Burkert et al. 2010), yielding $V_{\mathrm{circ}}(R_{E})=\sqrt{V(R_{E})^{2}+3.35\sigma_{V,0}^{2}}$ at $R_{E}$. We then calculate the total dynamical masses as

[TABLE]

where $G$ is the gravitational constant. We adopt a virial coefficient of $k_{\mathrm{tot}}(R_{E})=2.128$, which corresponds to a flattened Sérsic mass profile for a system with intrinsic axis ratio $q=0.4$ (to approximate the composite oblateness of the galaxy disk and halo; see Miller et al. 2011, Dutton et al. 2011a) and $n=1$ (exponential; typical for star-forming galaxies; Wuyts et al. 2011a), following the calculations given by Noordermeer (2008).444We explicitly apply an asymmetric drift correction and a constant virial coefficient $k_{\mathrm{tot}}(R_{E})$ for consistency with other work (e.g., Wuyts et al. 2016, Übler et al. 2017), as opposed to using $V_{\mathrm{rms}}(R_{E})$ and defining a composite $k_{\mathrm{eff}}(R_{E})$ as in Price et al. (2016). We note that this alternative approach would not impact the results of our previous study, with only a $\Delta\log_{10}M_{\mathrm{dyn}}=-0.006\mathrm{\,dex}$ mass difference between these methods. The impact of the assumed virial coefficient is discussed in Section 5.1.1.

For the resolved/aligned galaxies with detected rotation, we use the best-fit values of $V(R_{E})$ and $\sigma_{V,0}$ to directly determine $V_{\mathrm{circ}}(R_{E})$ and $M_{\mathrm{dyn}}$ using Equations 1 and 2.

3.3.2 Unresolved Kinematics from 1D Spectra

For galaxies without robustly detected rotation or that are unresolved/misaligned, we cannot simultaneously constrain $V(R_{E})$ and $\sigma_{V,0}$, as we measure their kinematics from integrated 1D spectra. We must instead assume a fixed value of $(V/\sigma_{V,0})_{R_{E}}$ to determine $V_{\mathrm{rms}}(R_{E})_{\mathrm{1D,corr}}$. From this, we can then calculate $V_{\mathrm{circ}}(R_{E})|_{(V/\sigma_{V,0})}$ (from Eq. 1):

[TABLE]

assuming an exponential disk profile, where $R_{E}=1.676\,r_{s}$. Finally, we calculate $M_{\mathrm{dyn}}$ for this fixed $(V/\sigma_{V,0})_{R_{E}}$ using Equation 2.

However, we   can   statistically   constrain   the   average $(V/\sigma_{V,0})_{R_{E}}$ for these galaxies by examining the correlation between the difference between dynamical and baryonic ($M_{\mathrm{baryon}}=M_{*}+M_{\mathrm{gas}}$) masses as a function of galaxy axis ratio. For disks, low inclination reduces the observed rotation signal while having no impact on an isotropic dispersion. For example, if $(V/\sigma_{V,0})_{R_{E}}$ is overestimated, $V(R_{E})$ will be overestimated and thus the composite $V_{\mathrm{rms}}(R_{E})$ will be overcorrected for inclination effects, leading to an overestimate of $M_{\mathrm{dyn}}$. For this case the degree to which $M_{\mathrm{dyn}}$ is overestimated should depend on axis ratio; edge-on ($b/a\sim 0.19$) galaxies have no inclination correction, while the correction is very high for face-on galaxies ($b/a\sim 1$). Since we expect no intrinsic trend between $M_{\mathrm{dyn}}$ or $M_{\mathrm{baryon}}$ and inclination (a random orientation projection effect), we would see this inclination overcorrection trend imprinted as a positive trend between $\Delta\log_{10}M=\log_{10}(M_{\mathrm{dyn}}/M_{\mathrm{baryon}})$ and $b/a$. Thus, we can infer the best-fit $V/\sigma_{V,0}$ as the value that best removes any dependence of $\Delta\log_{10}M$ on $b/a$. We note that this approach of modeling galaxies as inclined disks does not prevent us from finding very low to no ensemble-average rotation.

We constrain the ensemble average $(V/\sigma_{V,0})_{R_{E}}$ using the trend of $\Delta\log_{10}M$ versus $b/a$ as introduced in Price et al. (2016), but following a slightly different procedure. First, we calculate the dynamical masses for the galaxies without detected rotation and those that are unresolved/misaligned over a range of $(V/\sigma_{V,0})_{R_{E}}$ values (from $(V/\sigma_{V,0})_{R_{E}}=0$ to 10). We then measure the slope between $b/a$ and $\Delta\log_{10}M$ for these galaxies at every value of $(V/\sigma_{V,0})_{R_{E}}$, and determine the slope error by generating 500 bootstrap samples with replacement and measuring the slope for each realization. The best-fit $(V/\sigma_{V,0})_{R_{E}}$ is then estimated by finding the value that removes the trend of $\Delta\log_{10}M$ with $b/a$ (i.e., a slope of zero). The confidence interval is taken to be the range of $(V/\sigma_{V,0})_{R_{E}}$ where the slopes are consistent with zero within the uncertainties. Average $(V/\sigma_{V,0})_{R_{E}}$ for the unresolved/misaligned or undetected rotation galaxies are measured using this method in bins of redshift, stellar mass, sSFR, effective radius, and baryon surface density. The set of bin boundaries used in this work are shown with gray dashed lines in Figure 3.

We then use the best-fit $(V/\sigma_{V,0})_{R_{E}}$ in each bin to calculate $V_{\mathrm{rms}}(R_{E})_{\mathrm{1D,corr}}$ (Sec. 3.2), $V_{\mathrm{circ}}(R_{E})|_{(V/\sigma_{V,0})}$ (Eq. 3), and finally $M_{\mathrm{dyn}}$ (Eq. 2) for the nonresolved, non-rotation-detection galaxies. We additionally determine $(V/\sigma_{V,0})_{2.2}$ (see Sec 4.1.1) from the best-fit $(V/\sigma_{V,0})_{R_{E}}$, using the assumed arctan curve (Sec 3.2). The process is repeated for each of the binning parameter spaces ($z$, $z$-$M_{*}$, $z$-sSFR, $z$-$R_{E}$, and $z$-$\Sigma_{\mathrm{bar},e}$), allowing us to use consistently measured $(V/\sigma_{V,0})_{2.2}$ and $M_{\mathrm{dyn}}$ for the non-rotation-detected galaxies when examining trends in these parameter spaces.

3.3.3 Full-sample Dynamical Masses

The dynamical, baryonic, and stellar masses of the galaxies in our sample are shown in Figure 4. For the galaxies that are spatially unresolved/misaligned or have unresolved kinematics, the average $(V/\sigma_{V,0})_{R_{E}}$ is determined in bins of $z$ and $M_{*}$ (see Figure 3a). Overall, as expected, the dynamical masses exceed the stellar masses for most galaxies in all three redshift ranges, with a slightly higher offset at $z\sim 1.5$ relative to the higher redshift bins. The dynamical and baryonic masses are also reasonably consistent, with the median $M_{\mathrm{dyn}}$-$M_{\mathrm{baryon}}$ offset decreasing with increasing redshift. All galaxies (resolved/aligned galaxies with and without detected rotation, as well as the unresolved/misaligned galaxies) follow the same $M_{*}$-$M_{\mathrm{dyn}}$ and $M_{\mathrm{baryon}}$-$M_{\mathrm{dyn}}$ relations. The three subsamples also have similar distributions of $M_{\mathrm{dyn}}/M_{\mathrm{baryon}}$ versus axis ratio (i.e., inclination). We do note that our star-forming galaxy sample is incomplete below $\log_{10}(M_{\mathrm{baryon}}/M_{\odot})\lesssim 9.5-9.8$. The impact of this incompleteness is discussed further in Section 5.2.

4 Results

We use these kinematic and structural observations for galaxies in the MOSDEF survey to explore connections between kinematics and other properties and to constrain the evolution of the dynamical-baryonic mass offset and the inferred baryonic and dark matter fractions in galaxies between $z\sim 1.4-3.8$. Previous studies with more detailed kinematic observations have examined these properties, but the large sample size, three redshift intervals, and wide range of galaxy masses within each epoch of our sample provide a unique opportunity to investigate the kinematics and mass budgets for the star-forming galaxy population at $z\sim 1-3$.

4.1 Comparison between Internal Kinematics

and Galaxy Properties

We begin by examining how the kinematics of our sample vary with stellar mass ($M_{*}$), sSFR (which correlates with gas fraction; Tacconi et al. 2013), and baryon surface density ($\Sigma_{\mathrm{bar},e}$), in order to probe how galaxy structure and turbulence are set.

4.1.1 Trends with $V/\sigma_{V,0}$

The amount of kinematic support from ordered versus random motions ($V/\sigma_{V,0}$) provides information about the internal structures of galaxies. In particular, low values of $V/\sigma_{V,0}$ may indicate that a galaxy has a thick disk and high gas turbulence, while galaxies with high $V/\sigma_{V,0}$ tend to have ordered, thin disks. We thus investigate the relationship between $V/\sigma_{V,0}$ and other properties to constrain what processes drive the internal structures of star-forming galaxies at high redshifts.

In Figure 5, we present $(V/\sigma_{V,0})_{2.2}=V(2.2r_{s})/\sigma_{V,0}$555We consider $V(r=2.2r_{s})$, the radius where an exponential rotation curve peaks, for the analysis of the ratio of ordered-to-disordered motions ($V/\sigma_{V,0}$), to provide reasonable comparisons with existing measurements, as we do not directly constrain the turn-over or flattening for the rotation curves of our galaxies.

versus stellar mass (Figure 5a), sSFR (Figure 5b), baryon mass surface density, $\Sigma_{\mathrm{bar},e}$ (measured within the deprojected, major-axis effective radius $R_{E}$; Figure 5c), and redshift (Figure 5d) for the resolved/aligned galaxies. Median $(V/\sigma_{V,0})_{2.2}$ for the rotation-detected galaxies (large circles) are measured within bins of stellar mass, sSFR, and $\Sigma_{\mathrm{bar},e}$ split by redshift, and also binned by redshift alone, following the bin boundaries shown in Figure 3. The uncertainties on the medians are estimated by perturbing $(V/\sigma_{V,0})_{2.2}$ for each rotation-detected galaxy by its error over 500 realizations. We also show the ensemble average $(V/\sigma_{V,0})_{2.2}$ of the unresolved/misaligned and velocity-limit galaxies (diamonds). These ensemble-average values and their uncertainties are determined following the procedure of Section 3.3, and are separately measured within each of the above bins (again, see Figure 3). For clarity, the markers for best-fit median bins with no rotation ($(V/\sigma_{V,0})_{2.2}=0$) are displayed at $(V/\sigma_{V,0})_{2.2}=0.1$, given the log plot scale. Uncertainties that are consistent with the lower or upper ensemble fit boundary ($(V/\sigma_{V,0})_{2.2}=0$ or $10$) are marked with arrows. For reference, the median $(V/\sigma_{V,0})_{R_{E}}$ and $(V/\sigma_{V,0})_{2.2}$ for the resolved/aligned, rotation-detected galaxies and the unresolved/misaligned and rotation-limit galaxies within all bins are given in Tables 2 and 3, respectively.

We find a moderate positive correlation (Spearman rank correlation coefficient $\rho=0.30$ with $\sim\!\!2.9\sigma$ significance) between stellar mass and $(V/\sigma_{V,0})_{2.2}$ for our entire sample of galaxies with detected rotation. We observe a similar trend in the $(V/\sigma_{V,0})_{2.2}$ values measured in bins of $z$ and $M_{*}$ for the unresolved/misaligned and rotation-limit galaxies. Additionally, there is a moderate negative correlation ($\rho=-0.33$ with $>\!\!3\sigma$ significance) between $(V/\sigma_{V,0})_{2.2}$ and sSFR for the detected rotation galaxies, with a similar trend for the bins of unresolved/misaligned and rotation-limit galaxies. In contrast, we find a very weak, insignificant negative trend between $\Sigma_{\mathrm{bar},e}$ and $(V/\sigma_{V,0})_{2.2}$.

These trends of $(V/\sigma_{V,0})_{2.2}$ with $M_{*}$ and sSFR are in fairly good agreement with the findings of Newman et al. (2013) at $z\sim 1-2.5$, Wisnioski et al. (2015) at $z\sim 1$ and 2, Simons et al. (2016) and Alcorn et al. (2018) at $z\sim 2$, and Turner et al. (2017) at $z\sim 3.5$. We note that Wisnioski et al. find a larger trend of $V/\sigma_{V,0}$ with redshift at fixed stellar mass than we do, but the redshift intervals differ between the two studies ($z\sim 1$ and 2 versus $z\sim 1.5$, 2, and 3 used here). Additionally, the Wisnioski et al. sample excludes galaxies with $V/\sigma_{V,0}<1$, while ours does not.

Wisnioski et al. (2015) suggest that the trend of $V/\sigma_{V,0}$ with sSFR is indicative of galaxies tending to have more ordered motions when they have lower gas fractions, as predicted by the Toomre disk stability criterion (Toomre 1964; see also Genzel et al. 2011, Wisnioski et al. 2015). As sSFR correlates with gas fraction (Tacconi et al., 2013), we would thus expect to see a negative correlation between $V/\sigma_{V,0}$ and sSFR.

Since sSFR and stellar mass are also correlated, we would expect to observe both a decreasing trend of sSFR and an increasing trend of stellar mass toward higher $V/\sigma_{V,0}$. The slightly stronger correlation of $V/\sigma_{V,0}$ with sSFR that we observe suggests that gas fractions may be more important than stellar mass in setting a galaxy’s dynamical structure. Furthermore, the lack of correlation with $\Sigma_{\mathrm{bar},e}$ suggests that surface density has little impact on the relative amounts of ordered and unordered motions in these high-redshift galaxies.

The overall trends of $V/\sigma_{V,0}$ with $M_{*}$ and sSFR also appear to hold for galaxies at $z\sim 1.5$ and $z\sim 2$ separately. The trends with $V/\sigma_{V,0}$ at $z\sim 3$ are inconclusive due to the larger errors caused by small numbers of galaxies in some bins and the lack of dynamic range. Thus, while massive galaxies at $z\sim 1-2.5$ have thicker, more turbulent disks due to higher gas fractions than their local counterparts, the underlying physics regulating marginal disk stability may be the same and in place by $z\sim 2$.

For general reference, we also show the individual and median $(V/\sigma_{V,0})_{2.2}$ versus redshift (Figure 5d) for all galaxies, without controlling for any sample differences between the redshift bins. When we consider galaxies at fixed stellar mass, sSFR, or baryonic mass surface density, there are offsets that are suggestive of a decrease in $(V/\sigma_{V,0})_{2.2}$ with increasing redshift. However, given the measurement uncertainties, these offsets are also consistent with no evolution with redshift. The apparent decrease of $(V/\sigma_{V,0})_{2.2}$ with increasing redshift in Figure 5d thus reflects the average increase in sSFR (i.e., increase in gas fraction) and the average decrease in stellar mass between the redshift ranges (see Figure 1 and Tables 2, 3) and should not be interpreted as a pure redshift evolution. While our results show no conclusive redshift trends of $(V/\sigma_{V,0})_{2.2}$ with fixed galaxy mass, sSFR, or surface density, we cannot definitively rule out any $V/\sigma_{V,0}$ evolution over time.

Our median $V/\sigma_{V,0}$ values are in good agreement with the trend of previous results out to $z\sim 3.5$ (e.g., Epinat et al. 2008, Förster Schreiber et al. 2009, Law et al. 2009, Gnerucci et al. 2011, Epinat et al. 2012, Green et al. 2014, Wisnioski et al. 2015, Harrison et al. 2017, Swinbank et al. 2017, Turner et al. 2017; see Fig. 7 of Turner et al. 2017). As discussed above, our values for $(V/\sigma_{V,0})_{2.2}$ at $z\sim 2$ are somewhat lower compared to the KMOS${}^{\hbox{\scriptsize{3D}}}$ values because we do not apply a $V/\sigma_{V,0}$ cut. Additionally, our $z\sim 2$ sample extends to lower masses than the Wisnioski et al. $z\sim 2$ sample, which could further explain our lower $V/\sigma_{V,0}$ values.

4.1.2 Trends with Rotation and Intrinsic Dispersion for Galaxies with Robustly Detected Rotation

We also separately investigate correlations between rotation velocity, velocity dispersion, and other galaxy properties (as motivated in the beginning of Section 4.1) to directly examine how disk rotation velocity and turbulence evolve. Figure 6 shows measurements of $V_{2.2}=V(r=2.2r_{s})$ and $\sigma_{V,0}$ versus stellar mass, sSFR, baryon surface density, and redshift for the resolved/aligned galaxies with detected rotation. Median values within bins of redshift, stellar mass, sSFR, and $\Sigma_{\mathrm{bar},e}$ are also shown and are listed in Table 2.

We find the strongest correlation for our resolved/aligned rotation-detected sample between $V_{2.2}$ and $M_{*}$ ($\rho=0.48$ at $>\!\!3\sigma$). Additionally, there is fairly strong trend between $\sigma_{V,0}$ and $\Sigma_{\mathrm{bar},e}$ ($\rho=0.41$ at $>\!\!3\sigma$). There are also weaker, less significant correlations between $\sigma_{V,0}$ and $M_{*}$ ($\rho=0.29$ at $\sim\!2.8\sigma$) and $V_{2.2}$ and sSFR ($\rho=-0.28$ at $\sim\!2.7\sigma$). We do not find any correlation between $\sigma_{V,0}$ and sSFR or between $V_{2.2}$ and $\Sigma_{\mathrm{bar},e}$ for the rotation-detected sample.

These trends (or lack thereof) of $V_{2.2}$ and $\sigma_{V,0}$ are in qualitative agreement with previous work on massive star-forming galaxies at high redshifts. The observed correlation between $V_{2.2}$ and $M_{*}$ in the redshift slices agrees generally with the results presented by Harrison et al. (2017) at $z\sim 1$, Simons et al. (2016), Straatman et al. (2017), and Alcorn et al. (2018) at $z\sim 2$, and Turner et al. (2017) at $z\sim 3.5$. Similarly, the weak trend between $\sigma_{V,0}$ and $M_{*}$ split by redshift is consistent with the findings of Stott et al. (2016) at $z\sim 1$, Wisnioski et al. (2015) at $z\sim 1-2$, Alcorn et al. (2018) at $z\sim 2$, and Turner et al. (2017) at $z\sim 3.5$. In contrast, Simons et al. (2017) report no significant difference in dispersion with mass at redshifts between $z\sim 0.2-2$. However, we note that our median $\sigma_{V,0}$ values in bins of $M_{*}$ are consistent with the Simons et al. median values at both $z\sim 1.5$ and 2. Finally, the lack of correlation between $\sigma_{V,0}$ and sSFR that we see is similar to what Wisnioski et al. (2015) find.

The strong correspondence between rotation velocity and stellar mass is expected, as kinematics trace the total potential of a system. From the Toomre disk stability criterion, we would also expect the intrinsic dispersion to scale with both rotation velocity and gas fraction (e.g., Wisnioski et al. 2015). Additionally, the equilibrium or regulator galaxy growth model – where star-forming galaxies are in more or less steady equilibrium between gas inflow, star formation, and outflows (e.g., Bouché et al. 2010, Davé et al. 2012, Lilly et al. 2013, Dekel & Mandelker 2014) – would also predict a relation between intrinsic galaxy dispersion and gas fraction. Our finding of a slight correlation between $\sigma_{V,0}$ and stellar mass is in line with expectations from disk stability theory. However, we do not see a correlation between dispersion and sSFR (as a proxy for global gas fraction; e.g., Tacconi et al. 2013), as we would expect from the equilibrium model. Instead, we find that the intrinsic dispersion $\sigma_{V,0}$ is most strongly correlated with baryon surface density. The stronger correlation of dispersion with surface density, rather than sSFR or stellar mass, suggests that density may be more important in setting the amount of random motion than a galaxy’s global properties (e.g., total gas fraction or rotation velocity/total mass). However, given the uncertainties in our kinematic and sSFR measurements (and the inherent selection effects in our “rotation-detected” sample), our results cannot definitively rule out a relation between intrinsic velocity dispersion and gas fraction.

There is a general trend of $V_{2.2}$ increasing toward later times. This trend may primarily reflect the differences in stellar mass (and sSFR) in each redshift range (see Figure 3 and Table 2). When examining galaxies at fixed stellar mass, sSFR, and $\Sigma_{\mathrm{bar},e}$, we see do suggestive trends of increasing median $V_{2.2}$ over time. However, these offsets are also generally consistent with no redshift evolution, given the uncertainties and small number of “rotation-detected” galaxies. Nonetheless, considering the mass differences between redshift slices, our observed trends of $V_{2.2}$ with $z$ are generally consistent with the evolution from $z\sim 0.2$ to $z\sim 2$ found using DEEP2 and SIGMA data (Kassin et al. 2012, Simons et al. 2017) and agree with the overall trend of results from $z\sim 0$ out to $z\sim 3.5$ (Turner et al. 2017, and the literature values presented therein).

In comparison, we find suggestive trends of $\sigma_{V,0}$ decreasing from $z\sim 3.5$ to $z\sim 1.5$ when fixing stellar mass, sSFR, or $\Sigma_{\mathrm{bar},e}$, but these trends when holding galaxy properties fixed are also broadly consistent with no evolution given the uncertainties. Our observed values of $\sigma_{V,0}$ at $z\sim 1.4-3.3$ are in excellent agreement with the findings of DEEP2/SIGMA (Kassin et al. 2012, Simons et al. 2017) and are generally consistent with a range of studies at $z\sim 0-3.5$ (e.g., Newman et al. 2013; see Wisnioski et al. 2015, Turner et al. 2017 and included literature references). Even though we do not definitively detect a redshift trend in $\sigma_{V,0}$, the combination of our data and these previous measurements agrees fairly well with the expected evolution for marginally stable disks (e.g., Toomre 1964, Förster Schreiber et al. 2006, Genzel et al. 2011; see Fig. 8 of Wisnioski et al. 2015, Fig. 5 of Turner et al. 2017).

4.2 Evolution of the Dynamical-Baryonic Mass Offset

Next, by comparing the dynamical and baryonic masses for the galaxies in our sample, we investigate how the inferred baryonic and dark matter fractions vary with galaxy properties and evolve over time. Here we examine trends with stellar mass ($M_{*}$), effective radius ($R_{E}$), and baryon mass surface density ($\Sigma_{\mathrm{bar},e}$). These parameters enable us to examine how baryon and dark matter fractions vary over redshift as a function of galaxy mass and size (as a proxy for the extent into the halo). We will also consider whether there is a redshift-independent relation between baryon surface density and total baryon fraction (as suggested by Wuyts et al. 2016).

The sample is split into bins of redshift with stellar mass, effective radius, and baryon mass surface density as shown in Figure 3. The ensemble average $(V/\sigma_{V,0})_{R_{E}}$ for the galaxies without detected rotation and that are unresolved/misaligned are measured within these bins following the procedure described in Section 3.3. We then combine all galaxies within each bin (all of the detected rotation, rotation-limit, and unresolved/misaligned galaxies) and calculate the median offset between dynamical and baryonic masses, $\Delta\log_{10}M=\log_{10}(M_{\mathrm{dyn}}/M_{\mathrm{baryon}})$, and determine errors through bootstrap resampling.

These mass offsets are shown versus redshift in Figure 7, where each bin is colored by the median $M_{*}$, $R_{E}$, and $\Sigma_{\mathrm{bar},e}$ (left to right, respectively). We find that the mass difference $\Delta\log_{10}M$ decreases toward higher redshifts, where galaxies have higher baryon fractions ($f_{\mathrm{bar}}=M_{\mathrm{baryon}}/M_{\mathrm{dyn}}$). This overall trend of increasing mass offset with time holds even in fixed bins of stellar mass, effective radius, and surface density, suggesting that this general result is not solely driven by differences in these properties for our sample within different redshift ranges. The offsets within these bins are also given in Table 4. These results agree well with previous work finding that massive galaxies are typically baryon dominated within the galaxy scale at high redshift (e.g., Förster Schreiber et al. 2009, Wuyts et al. 2016, Burkert et al. 2016, Stott et al. 2016, Alcorn et al. 2016, Genzel et al. 2017, Lang et al. 2017), though these studies did not probe the mass offsets as a joint function of redshift and other galaxy properties. Other studies of high-redshift galaxies find much higher dark matter fractions within much larger radii (e.g., $\gtrsim 60\%$ within $6r_{s}$; Tiley et al. 2019). While we do not have the S/N in the outskirts of our galaxies to constrain the dark matter fractions within such large radii (here we only constrain the values within $R_{E}$), these large-radius results are not necessarily inconsistent with the smaller-radius results, as the dark matter fraction should naturally be higher when probing apertures farther out into the halo.

This observed baryon fraction trend agrees fairly well with theoretical work. Results from hydrodynamical simulations also show an average increase in the baryon fraction of star-forming or disk galaxies within their half-light (or half-mass) radii toward higher redshifts (e.g., Lovell et al. 2018; $z\sim 0-4$, Teklu et al. 2018; $z\sim 2$). However, theory and observations begin to have tension at $z\gtrsim 2$, where the observed very high or unphysical $f_{\mathrm{bar}}$ are higher than the predicted baryon fractions from simulations (Wuyts et al. 2016, Lovell et al. 2018, Teklu et al. 2018).

We find little to no variation in mass offset with $M_{*}$ at fixed redshift (Figure 7a). At $z\sim 1.5$, the median mass offset of all $M_{*}$ bins corresponds to baryon fractions within the effective radius of $f_{\mathrm{bar}}\sim\!65\%$, or dark matter fractions of $f_{\mathrm{DM}}\sim\!35\%$ (where $f_{\mathrm{DM}}=1-f_{\mathrm{bar}}$). For all stellar mass bins at $z\gtrsim 2$, the mass offsets are negative and fall within the unphysical region where $M_{\mathrm{dyn}}<M_{\mathrm{baryon}}$. However, we note that the highest stellar mass bin ($\log_{10}(M_{*}/M_{\odot})>10.25$) at $z\sim 2.3$ is consistent with up to $f_{\mathrm{DM}}\sim\!5\%$ (or $f_{\mathrm{bar}}\sim\!95\%$) within the uncertainties. Our results at $z\sim 1.5$ and $z\sim 2.3$ for the highest stellar mass bin are in relatively good agreement with the results by Wuyts et al. (2016) for similar redshifts and masses. The relative lack of variation of $f_{\mathrm{bar}}$ with stellar mass observed in this work is in contrast to predictions from theoretical work. Results from the Illustris TNG simulations at $z\sim 2$ predict an average increase in the galaxy-scale (e.g., $r_{\mathrm{half}}$) baryon fraction for disk galaxies with increasing stellar mass (Lovell et al., 2018). Furthermore, Lovell et al. (2018) predict very little baryon or dark matter fraction evolution in galaxies with $\log_{10}(M_{*}/M_{\odot})=9$ from $z\sim 4$ to the present day, in contrast to the results we find for our lowest mass bin ($\log_{10}(M_{*}/M_{\odot})\sim 9.5$).

When controlling for galaxy size, we find larger galaxies have higher $\Delta\log_{10}M$ and lower baryon fractions (Figure 7b). For galaxies at $z\sim 1.5$, galaxies with $R_{E}\geq 2.5\mathrm{\,kpc}$ have $f_{\mathrm{bar}}\sim 55\%$ within their effective radii, while smaller galaxies with $R_{E}<2.5\mathrm{\,kpc}$ have $f_{\mathrm{bar}}>100\%$ and no inferred dark matter (but are consistent with up to a $\sim\!10\%$ dark matter fraction within the uncertainties). At $z\sim 2.3$ and $z\sim 3.3$, the baryon fractions for galaxies with $R_{E}\geq 2.5\mathrm{\,kpc}$ are unphysical ($f_{\mathrm{bar}}\gtrsim 100\%$) but are consistent with small dark matter fractions, while the smaller galaxies fall fully within the unphysical regime.

Finally, we examine the mass offset as a function of baryon surface density at fixed redshift (Figure 7c). We find that denser galaxies have lower $\Delta\log_{10}M$ and higher inferred baryon fractions. At $z\sim 1.5$, the lowest- and medium-density galaxies ($\langle\log_{10}(\Sigma_{\mathrm{bar},e}/[M_{\odot}\mathrm{\,kpc}^{-2}])\rangle\sim 8.4,8.8$) have $f_{\mathrm{bar}}\sim 50\%,75\%$ within $R_{E}$, respectively, while the highest-density galaxies ($\langle\log_{10}(\Sigma_{\mathrm{bar},e}/[M_{\odot}\mathrm{\,kpc}^{-2}])\rangle\sim 9.2$) have $f_{\mathrm{bar}}>100\%$, but are consistent with $f_{\mathrm{DM}}\sim 15-25\%$ within the uncertainties. The mass offsets decrease to higher redshifts, with the lowest-density galaxies having increased $f_{\mathrm{bar}}\sim 90\%$ at $z\sim 2.3$ and $f_{\mathrm{bar}}>100\%$ (no dark matter) at $z\sim 3.3$, and with he medium- and high-density galaxy bins all having $f_{\mathrm{bar}}>100\%$ at $z\gtrsim 2$.

While our results suggest that $f_{\mathrm{bar}}$ depends little on mass, the extent of a galaxy within its halo could explain variations of baryon fraction with size, density, and redshift (Price et al. 2016, Wuyts et al. 2016, Übler et al. 2017). Larger or less dense galaxies may extend farther into their halos, resulting in lower $f_{\mathrm{bar}}$ within their effective radii. The average size growth of galaxies over time would likewise be expected to lead to lower baryon densities at lower redshifts.

Our observed trends of increasing $f_{\mathrm{bar}}$ with increasing density ($\Sigma_{\mathrm{bar},e}$), decreasing radius ($R_{E}$), and increasing redshift support this postulate that galaxy extent within its halo impacts the observed $f_{\mathrm{bar}}$ within $R_{E}$. This finding is in excellent agreement with the full-sample results of Wuyts et al. (2016) ($z\sim 0.6-2.6$) and generally agrees with the correlation between $f_{\mathrm{bar}}$ and $\Sigma_{\mathrm{bar},e}$ seen for galaxies from the Illustris simulation matched to the KMOS${}^{\hbox{\scriptsize{3D}}}$ sample (Wuyts et al., 2016). However, our results also show that there is redshift evolution in the baryon fraction-density ($f_{\mathrm{bar}}$-$\Sigma_{\mathrm{bar},e}$) and -size ($f_{\mathrm{bar}}$-$R_{E}$) trends, suggesting that neither size nor density alone sets a galaxy’s baryon fraction.

A key factor we have not included is how the dark matter halo mass profiles evolve and how galaxy and halo properties are related. Exact theoretical predictions of central halo profiles and their evolution are uncertain, given the potential impact of adiabatic halo contraction (or expansion) and the response of the halo to a baryonic disk (e.g., see discussion in Duffy et al. 2010, Dutton & Macciò 2014, Velliscig et al. 2014, Courteau & Dutton 2015, Dutton et al. 2016). The particular details of how the halo and galaxy coevolve could strongly impact the interpretation of our derived central baryon/dark matter fraction results. For instance, if $R_{E}/R_{\mathrm{halo}}$ changes with redshift, holding $R_{E}$ fixed would also result in $\Delta\log_{10}M$ probing halo evolution within this fixed physical radius (e.g., concentration). Likewise, evolution in $M_{*}/M_{\mathrm{halo}}$, as well as halo mass profile evolution, could also impact the redshift trends of $\Delta\log_{10}M$ at fixed $M_{*}$ and $\Sigma_{\mathrm{bar},e}$. Taken together, these results suggest that the evolution of galaxies’ baryon fractions over time reflects the complex interplay of galaxy mass-size growth, global gas fraction, and halo growth and evolution, and thus it is not surprising that none of the individual examined properties are responsible for a universal, redshift-independent $f_{\mathrm{bar}}$ relation.

The negative offsets $\Delta\log_{10}M$ (i.e., $f_{\mathrm{bar}}\!>\!100\%$) observed at $z\gtrsim 2$, especially for smaller and denser galaxies, imply that at least some of the assumptions used to measure the dynamical and baryonic masses may be invalid for these galaxies. At minimum, we would expect the dynamical mass to account for all of the observed baryonic mass. This tension, as well as potential ways the masses might be reconciled, is further discussed in Section 5.1.1. Nonetheless, the observed redshift evolution of $f_{\mathrm{bar}}$ is likely not entirely driven by these mass measurement uncertainties. We would expect potential underestimates of dynamical masses or overestimates of gas masses to be the least problematic for larger, less dense galaxies (as they mostly do not have unphysical baryon fractions), and we do see time evolution of $f_{\mathrm{bar}}$ for these galaxies. In contrast, the exact details of the evolution of $f_{\mathrm{bar}}$/$f_{\mathrm{DM}}$ for smaller, denser galaxies will require more detailed future studies.

5 Discussion

5.1 Reconciling Baryonic and Dynamical Masses

at High Redshifts

While our results support a decreasing dark matter fraction within $R_{E}$ toward higher redshifts, there is tension between the dynamical and baryonic masses for a large fraction ($54\%$) of our sample, particularly in the higher redshift bins. On average, galaxies at $z\sim 1.5$ and for the highest mass bin at $z\sim 2$ are consistent with nonzero dark matter fractions. The $z\sim 1.5$ median mass offsets leave room for a more bottom-heavy IMF (e.g., Salpeter 1955), but other studies suggest that a Chabrier (2003) IMF is more appropriate for star-forming galaxies (e.g., Bell & de Jong 2001, Tacconi et al. 2008, Dutton et al. 2011b, Brewer et al. 2012). However, even for a Chabrier (2003) IMF, the dynamical masses on average are lower than the inferred baryonic masses for all stellar masses at $z\gtrsim 2$. Similarly, when splitting the sample by $R_{E}$ or $\Sigma_{\mathrm{bar},e}$, small and dense galaxies generally have lower inferred dynamical than baryonic masses. The unphysical offsets where $M_{\mathrm{dyn}}<M_{\mathrm{baryon}}$ suggest that some of the assumptions used to derive masses are invalid for these galaxies. We discuss and explore possible solutions to these tensions below.

5.1.1 Impact of Higher Star Formation Efficiencies or Virial Coefficients

Two potential causes of the tension between the dynamical and baryonic masses (or possibly both) are systematic overestimates of the baryonic masses and systematic underestimates of the dynamical masses. Given that we do not directly measure molecular gas masses, nor do we have the detailed, high-S/N observations required to fully model the dynamical masses of our galaxies, we have adopted prescriptions to infer the molecular gas masses from measured SFRs and to convert measured velocities and dispersions into dynamical masses. Both of our adopted prescriptions could therefore contribute to our observed mass tension.

To start, if we overestimate the galaxies’ gas masses, we would also overestimate their baryonic masses. In this study, we have used the SFR-gas mass relation of Kennicutt (1998) for star-forming galaxies in the local universe to convert observed, dust-corrected SFRs into gas masses. However, this relation may not hold for galaxies at higher redshifts, particularly for higher surface density galaxies. Work on galaxies at $z\sim 1-3$ finds Kennicutt-Schmidt relation slopes that vary from $N=1.28$ (Genzel et al., 2010) to $N=1.7$ (Bouché et al., 2007), which bracket the local value of $N=1.4$ by Kennicutt (1998). If the true slope for our galaxies is higher than the local relation (i.e., closer to $N=1.7$), our inferred gas masses would overestimate the true values, which would ease the tension between the dynamical and baryonic masses.

In Figure 8a, we explore the possibility that the some gas masses may be overestimated by examining what gas surface densities (i.e., modified star formation efficiencies [SFEs]) would be needed to have consistent masses (i.e., $M_{\mathrm{dyn}}\geq M_{\mathrm{baryon}}$). We show the relation of Kennicutt (1998) with a thick black line, which we use to infer gas masses (through $\Sigma_{\mathrm{gas},e}$) from the SFR surface densities of our sample. We then calculate the maximum gas surface density $\Sigma_{\mathrm{gas},e}^{\mathrm{new}}$ the “unphysical” ($M_{\mathrm{dyn}}<M_{\mathrm{baryon}}$) galaxies could have and still have physically consistent masses, using $M_{\mathrm{gas}}^{\mathrm{new}}=M_{\mathrm{dyn}}-M_{*}$ (i.e., no dark matter). With the Kennicutt (1998)-derived gas masses, 54% of our sample has unphysical masses, with $M_{\mathrm{dyn}}<M_{\mathrm{baryon}}$. When calculating the maximum physically allowable gas masses for these unphysical baryon fraction galaxies (i.e., higher SFE and decreased gas surface densities), we find that 76% of our sample (up from 46% for our fiducial calculation) falls between the “normal” star-forming galaxy and ULIRG/SMG galaxy star-formation surface density to gas surface density regimes by Genzel et al. (2010) (marked with black and red dashed lines, respectively). Such an increase in the SFE for these galaxies may be reasonable, given observations showing a decrease in gas depletion time (i.e., increased SFE) toward higher redshifts (e.g., Tacconi et al. 2018). Even more extreme SFEs than the ULIRG/SMG relation by Genzel et al. (2010) would be needed to bring the masses into agreement for the remaining galaxies with $M_{\mathrm{dyn}}\geq M_{*}$ (10% of total sample). SFE changes cannot be invoked for the galaxies with $M_{*}>M_{\mathrm{dyn}}$ (14% of total sample), though given the uncertainties, we may expect some scatter even into this region. Overall, in agreement with a similar test by Wuyts et al. (2016), this check suggests that moderately increased SFEs relative to the Kennicutt (1998) Schmidt relation could explain the mass tension for the majority of our sample, but further work and observations are needed to fully quantify the Schmidt relation in high-redshift galaxies.

Alternatively (or additionally), our measurements and prescriptions may underestimate the true dynamical masses of our sample. First, if the ionized gas kinematics do not fully trace the total system potential, then the dynamical masses inferred from the kinematics will underestimate the true total masses. Second, we have assumed an underlying oblate intrinsic mass profile and a constant mass-to-light ratio in order to infer total dynamical masses from the measured kinematics at $R_{E}$, through the virial coefficient $k_{\mathrm{tot}}$. If a galaxy’s true mass profile is more spherical than the assumed oblate mass profile, using the adopted $k_{\mathrm{tot}}$ can underestimate the total galaxy mass.

Similar to the test with gas masses, we explore what virial coefficient $k_{\mathrm{tot}}$ would need to be adopted to reconcile the masses in Figure 8b. We show the fiducial $k_{\mathrm{tot}}=2.128$ as a thick black horizontal line. For objects with $M_{\mathrm{dyn}}<M_{\mathrm{baryon}}$, we calculate the minimum modified $k_{\mathrm{tot}}$ that would be necessary for physical masses, using the limiting case of $M_{\mathrm{dyn}}=M_{\mathrm{baryon}}$, or taking $k_{\mathrm{new}}=M_{\mathrm{baryon}}G/(V_{\mathrm{circ}}^{2}R_{E})$. After determining $k_{\mathrm{new}}$ for the unphysical baryon fraction galaxies, 62% of the galaxies (up from 46%) would have virial coefficients falling between the fiducial $k_{\mathrm{tot}}$ for a $q=0.4$ oblate exponential profile and the value for a spherical exponential potential ($k_{\mathrm{tot}}=2.933$, dotted black horizontal line). The remainder would need even higher virial coefficients (corresponding to prolate potentials), up to very extreme, unrealistic values ($k_{\mathrm{tot}}\gtrsim 10$). van der Wel et al. (2014b) have argued that a large fraction of massive galaxies at $z\gtrsim 2$ have spheroidal/disk geometries, while the majority of low-mass galaxies at $z\gtrsim 2$ are elongated and prolate, using an analysis of the axis ratio distribution of galaxies out to $z\sim 2.5$, which might suggest virial coefficients higher than the oblate value we have used in this work. We do note that very few galaxies in our kinematics sample at $z\sim 3$ have axis ratios $b/a>0.8$, as expected for prolate geometries, but this may only reflect the relatively small sample size at $z\sim 3$. Nonetheless, much more detailed kinematic observations of small galaxies at $z\sim 2-3$ would be necessary to fully constrain their internal geometries.

To further quantify the impact of a higher SFE or virial coefficient, we also determine how the median dynamical-to-baryonic mass ratio, $\Delta\log_{10}M$, changes if we adopted the ULIRG/SMG relation from Genzel et al. (2010) or the spherical $k_{\mathrm{tot}}$ for our full sample. We find that these modifications would increase $\Delta\log_{10}M$ by $\sim\!0.27$ and $\sim\!0.14\mathrm{\,dex}$, respectively. Figure 9 shows these offsets, along with other systematic uncertainties discussed throughout Section 5.

5.1.2 Alternative Gas Mass Estimates: Scaling Relations

We next examine how our mass fraction results (e.g., $f_{\mathrm{bar}}$) would change — and whether any of the dynamical-baryonic mass tension could be alleviated — if we measure gas masses for our full sample following high-redshift gas mass scaling relations instead of using the Kennicutt (1998) Schmidt relation. For this test, we estimate gas masses using the best-fit molecular gas scaling relation of Tacconi et al. (2018), using the Whitaker et al. (2014) main-sequence prescription. We also adopt the “best” SFRs (and matching sSFRs) derived following the ladder technique of Wuyts et al. (2011b) in place of our primarily $\textrm{H}\alpha$- and $\textrm{H}\beta$-based SFRs, as these SFRs match those used to derive the scaling relations. We then remeasure $(V/\sigma_{V,0})_{R_{E}}$ for the galaxies without detected rotation using the scaling-relation baryonic masses (following Section 3.3.2), and redetermine the dynamical-baryonic mass offset as in Section 4.2.

Our primary results are unchanged if scaling-relation gas masses are used. On average the scaling-relation gas masses are slightly higher $\sim 0.13\mathrm{\,dex}$) than the Schmidt-relation gas masses, resulting in even higher average $f_{\mathrm{bar}}$ (the median $\Delta\log_{10}M$ decreases by $\sim\!-0.07\mathrm{\,dex}$; see Fig. 9). Nonetheless, the general trend of increasing $f_{\mathrm{bar}}$ toward higher redshifts at fixed mass, size, or surface density is unchanged, and the dynamical-baryonic mass tension at $z\gtrsim 2$ is still present when adopting these alternative gas masses. When using the scaling-relation gas masses, there is a suggestive trend in our lowest redshift bin ($z\sim 1.5$) where higher stellar mass galaxies have slightly higher $f_{\mathrm{bar}}$. However, the current data uncertainties hamper a firm conclusion regarding any definitive differences between this finding and our fiducial results where $f_{\mathrm{bar}}$ does not depend on $M_{*}$ at a given redshift. Overall, we conclude that our results are not greatly impacted by the choice of adopting gas masses estimated by inverting the Kennicutt (1998) relation instead of masses derived using scaling relations.

5.1.3 Impact of Incompleteness at Low Baryonic Masses

Another concern regarding our dynamical-baryonic mass comparison comes from the incompleteness of our sample at the lowest baryonic masses ($\log_{10}(M_{\mathrm{baryon}}/M_{\odot})\lesssim 9.5-9.8$; see Figure 4). This incompleteness can lead to underestimates of $\Delta\log_{10}M$, particularly at the lowest masses, as objects with lower $M_{\mathrm{baryon}}$ at fixed $M_{\mathrm{dyn}}$ near the completeness limit are missing from this analysis. This underestimate of $\Delta\log_{10}M$ is equivalent to an overestimate of $f_{\mathrm{bar}}$, especially at lower stellar masses and higher redshifts, as the median stellar mass of our sample decreases toward higher redshifts. We examine the impact of this incompleteness by repeating the analysis of Section 4.2, but including only galaxies with $\log_{10}(M_{\mathrm{dyn}}/M_{\odot})\geq 9.8$. We do not refit the unresolved/misaligned/rotation nondetection galaxies $(V/\sigma_{V,0})_{R_{E}}$, but simply apply the cut to the fiducial $M_{\mathrm{dyn}}$.

We find increases of $\sim\!0.17$, $\sim\!0.19$, and $\sim\!0.15\mathrm{\,dex}$ in $\Delta\log_{10}M$ for the lowest stellar mass bin ($\log_{10}(M_{*}/M_{\odot})<9.75$) at $z\sim 1.5$, $z\sim 2$, and $z\sim 3$, respectively, when applying this dynamical mass cut. Altogether, the full-sample median mass offset increases by $\sim\!0.08\mathrm{\,dex}$ (see Fig. 9). However, the greater impact of incompleteness at lower masses could be driving part of the tension with theoretical results. Based on the shifts from this dynamical mass cut, accounting for incompleteness could result in lower baryon fractions for the lower-mass galaxies relative to the higher-mass galaxies, as found by Lovell et al. (2018). Nonetheless, our primary results remain unchanged: $f_{\mathrm{bar}}$ increases toward higher redshifts, and most galaxies still have $M_{\mathrm{dyn}}<M_{\mathrm{baryon}}$ at $z\gtrsim 2$. We conclude that, on the whole, incompleteness does not solely drive our results but could be partially responsible for the lack of offset in $f_{\mathrm{bar}}$ between stellar mass bins.

5.1.4 Unresolved Kinematics: Early-type Galaxies?

For this analysis, we have modeled the galaxies that are unresolved, are misaligned, or have no detected rotation as inclined thick disks, but this may not hold for all of these objects. Instead, it is possible that the dynamics of at least some of these systems are better described by early-type, random-motion-dominated structures. While our models and adopted dynamical mass calculation do allow for no ensemble rotation ($(V/\sigma_{V,0})_{R_{E}}=0$), resulting in inclination-independent values and an “effective” virial coefficient of $k_{\mathrm{tot}}=7.13$, there may still be differences in our fiducial $M_{\mathrm{dyn}}$ and masses determined following prescriptions for early-type galaxies. As we cannot directly distinguish between inclined thick disks and early-type dynamics in our integrated 1D velocity dispersion measurements, we thus also calculate the dynamical masses of all the unresolved/misaligned and rotation-limit galaxies under the assumption that they are early-type galaxies (following Cappellari et al. 2006; see Section 5.2 of Price et al. 2016), to examine how the assumed kinematic structure of these galaxies impacts our results.

To determine early-type dynamical masses, we calculate aperture-corrected integrated dispersions $\sigma_{e,\mathrm{corr}}$ from the observed, instrument-corrected 1D velocity dispersions $\sigma_{V,\,\mathrm{obs}}$. We then use circularized radii to calculate both the dynamical and gas masses for consistency. As we found in our previous study, the dynamical and baryonic masses for the full sample are still in reasonable agreement ($\Delta\log_{10}M=0.06,-0.18,-0.23$ at $z\sim 1.5,2,3$, respectively; the full-sample median offset is $\sim\!-0.07\mathrm{\,dex}$, see Fig. 9), and there is no strong residual trend of $\Delta\log_{10}M$ with $b/a$. However, as we found before, under the early-type assumption, the unresolved/misaligned/ rotation-undetected galaxies have systematically lower

$\Delta\log_{10}M$ than the rotation-detected galaxies, with even higher inferred baryon fractions. This increased dynamical-baryonic mass tension suggests that, on average, these galaxies have at least a modest amount of rotational support, supporting the adoption of the best-fit ensemble average $(V/\sigma_{V,0})_{R_{E}}$ and the corresponding $M_{\mathrm{dyn}}$ as derived in Section 3.3.2 for this study. However, systematic uncertainties (see Section 5.2) make it difficult to definitively conclude whether or not these galaxies have rotational support.

5.2 Other Analysis Caveats and the Path Forward

In addition to the caveats discussed in the previous section, other assumptions made in this work could impact our results. First, many of the caveats discussed in Price et al. (2016) also apply to this analysis. Specifically, the accuracy of the Galfit structural parameters is not fully accounted for in our kinematic measurements. These structural parameters are integral to both the 1D and 2D kinematic modeling, but no structural parameter errors are included in the kinematic fitting. Our analysis also will have over- and underestimates in the measured kinematics if the photometric and kinematic major axes are misaligned, resulting in different position angle misalignments $\Delta\mathrm{PA}$ for the slit and the photometric and kinematic major axes and introducing extra scatter into our measurements. Similarly, the uncertainty in the measured galaxy axis ratios impacts the inclination correction, particularly for smaller galaxies (e.g., $z\sim 3$ galaxies). If the observed ratios are larger or smaller than the true values, then this will lead to over- and underestimates in the kinematics and dynamical masses, respectively, which will also increase the scatter of our measurements. Additionally, we assume a fixed intrinsic disk thickness of $(b/a)_{0}=0.19$ to derive galaxy inclinations. Variations in the true intrinsic thicknesses (e.g., thicker or thinner) will lead to over- and underestimates (respectively) of the inclination correction for our galaxies, which could increase the scatter and may also introduce systematic offsets in the derived dynamical properties. For thicker $(b/a)_{0}=0.3$ or thinner $(b/a)_{0}=0.1$, the derived dynamical masses could be up to $\sim-0.025\mathrm{\,dex}$ lower and $\sim 0.011\mathrm{\,dex}$ higher than the fiducial dynamical masses, respectively (see Fig. 9).

Furthermore, we have assumed fixed forms for the kinematic profiles of our galaxies. We adopt an arctan rotation curve model (e.g., following Weiner et al. 2006, Miller et al. 2011). Some work suggests that particularly massive high-redshift star-forming galaxies have falling rotation curves, such as a Freeman (1970) exponential disk model (e.g., Wisnioski et al. 2015, van Dokkum et al. 2015, Genzel et al. 2017, Lang et al. 2017). However, for our objects, we only reliably probe the kinematics out to $\sim 2.2r_{s}$. Arctan model profiles are similar to the other rotation curve models over this range, so our choice of rotation profile should not impact our results very strongly. We also assume a constant intrinsic velocity dispersion profile. If the true velocity dispersion instead decreases with increasing radius, this would result in higher median $(V/\sigma_{V,0})_{R_{E}}$ measured for the galaxies without resolved rotation (also see discussion in Section 5.7 of Price et al. 2016). Higher spatial resolution and higher-S/N observations are needed to properly constrain the form of the velocity dispersion profiles.

Second, beyond the assumption of a SFR-gas mass relation (e.g., from Kennicutt 1998, as discussed in Section 5.1.1), we also use a mix of SFR indicators, as Balmer-decrement-corrected $\textrm{H}\alpha$ SFRs are not available for the entire sample. Mixing SFR indicators could lead to systematic differences in the inferred gas masses. However, in Appendix B, we compare SFR indicators for galaxies at $z\sim 2$ with Balmer-decrement-corrected $\textrm{H}\alpha$ SFRs and find no bias between the $\textrm{H}\alpha$/$\textrm{H}\beta$ +$A_{\mathrm{V}}$ indicators and $\textrm{H}\alpha$ SFRs, though the SED SFRs underestimate the $\textrm{H}\alpha$ SFRs (as also found by Reddy et al. 2015 for galaxies with high SFRs). This offset would lead to even larger inferred gas masses, increasing the tension between the dynamical and baryonic mass measurements. If we use the median offset between the other SFR indicators and the $\textrm{H}\alpha$ SFRs (as measured in Appendix B) to correct all non-Balmer-decrement $\textrm{H}\alpha$ SFRs, we find the same median dynamical-to-baryonic mass ratio for the full sample (or at low or high masses), because the SED SFRs are only used for $\sim\!11\%$ of the sample.

However, some $z\sim 3$ galaxies with very high sSFRs may have the opposite problem, where the SED and $\textrm{H}\beta$ +$A_{\mathrm{V}}$ SFR indicators may overestimate the intrinsic $\textrm{H}\alpha$ SFRs (as discussed in Appendix B). As directly measured Balmer-decrement $\textrm{H}\alpha$ SFRs for $z\sim 3$ galaxies are currently unavailable, we instead test how $\Delta\log_{10}M$ changes at $z\sim 3$ if we only include $z\sim 3$ galaxies with $\log_{10}(\mathrm{sSFR}/\mathrm{\,yr}^{-1})\leq-7.8$. When performing this cut, we find a $\sim 0.07\mathrm{\,dex}$ increase in the mass offset at $z\sim 3$ (see Fig. 9), which is insufficient to move the median mass offsets at $z\sim 3$ into the physical regime.

Third, we have assumed that the galaxies have constant mass-to-light profiles (i.e., equal half-light and half-mass radii). If half-mass radii are smaller than rest-frame optical $R_{E}$, as found by Szomoru et al. (2013), we would overestimate the true dynamical masses, which would result in increased tension between the dynamical and baryonic masses. Specifically, for the average ratio $r_{1/2,\mathrm{mass}}/R_{E}\sim 0.75$ found by Szomoru et al. (2013), dynamical masses measured using half-mass radii would be lower than the fiducial $M_{\mathrm{dyn}}$ by $\sim\!-0.15\mathrm{\,dex}$ (see Fig. 9). Recently, Suess et al. (2019) found that the median ratio $R_{E,\mathrm{mass}}/R_{E,\mathrm{light}}$ for star-forming galaxies decreases from $\sim 1$ at $z\gtrsim 2$ down to $\sim 0.75-0.8$ at $z\sim 1.5$, so this reduction in $M_{\mathrm{dyn}}$ would likely impact our lowest redshift bin more than our higher redshift bins. Suess et al. also find that $R_{E,\mathrm{mass}}/R_{E,\mathrm{light}}$ decreases with stellar mass, which would imply that at high $M_{*}$ we would find lower $\Delta\log_{10}M$ (higher $f_{\mathrm{bar}}$), with less offset for lower masses. To quantify the impact variable mass-to-light profiles may have on our results, we use the redshift- and mass-dependent ratios between $R_{E,\mathrm{mass}}/R_{E,\mathrm{light}}$ from Suess et al. (2019) to calculate how $M_{\mathrm{dyn}}$ changes. We assume $r_{t}=0.4r_{s}=0.4R_{E}/1.676$ and use either the measured $V(R_{E})$ or inferred $V(R_{E})$ (from $(V/\sigma_{V,0})_{R_{E}}$ fit in bins of redshift and stellar mass), and we also apply the extrapolated relations to stellar masses below the Suess et al. sample completeness limits. We find essentially no change in the full-sample median dynamical-to-baryonic mass ratio (formally a $\sim\!0.004\mathrm{\,dex}$ increase; see Fig. 9), though for the mass and redshift bins where Suess et al. find ratios of $R_{E,\mathrm{mass}}/R_{E,\mathrm{light}}\sim 0.75-0.8$, the offset within those bins will be close to the $\sim\!-0.15\mathrm{\,dex}$ shift for the Szomoru et al. case.

This problem of smaller half-mass than half-light radii would be further compounded if emission-line half-light radii were adopted, since the emission line profiles generally extend farther out than the stellar light profiles (Nelson et al., 2016). In this case, the half-emission-line radii would be even larger than the half-mass radii. Since our dynamical mass derivation requires kinematics measured at the half-mass radius, if kinematics and sizes were instead measured at the on-average larger $R_{E,\textrm{H}\alpha}$, we would overestimate the dynamical masses of our sample.

However, in Price et al. (2016) we found that the seeing-matched $\textrm{H}\alpha$ and stellar light profiles for the resolved rotation galaxies were similar, which suggests that using stellar light profiles to model and measure galaxy kinematics is a reasonable assumption for most of our sample. In contrast, if the half-mass radii are larger than $R_{E}$ for some galaxies, the presented values would be underestimates of the total dynamical masses. Mock observations of multiple lines of sight to the same galaxy show that the specific viewing direction to a galaxy can strongly impact the observed half-light radius, leading to additional variations between the half-light and half-mass radii (including cases with half-light radii smaller than the half-mass radius; Price et al. 2017). This potential line-of-sight size variation would introduce additional scatter into the dynamical mass measurements.

Finally, as mentioned in Section 5.1.1, if the ionized gas kinematics do not fully trace the potential of these galaxies (at least out to the half-mass radius), the kinematics and dynamical masses will underestimate the true values. This concern would particularly impact lower-S/N objects and small galaxies close to the seeing limit, where we determine the kinematics from integrated 1D dispersion measurements.

Despite these analysis caveats, the following conclusions are fairly robust. First, while the precise baryon and dark matter fractions are somewhat uncertain, there is strong evidence that massive star-forming galaxies are highly baryon dominated at high redshifts. Furthermore, small and compact galaxies at high redshifts are more baryon dominated than more extended galaxies, suggesting that the extent of a galaxy within its halo plays a key role in setting the dark matter fraction within the half-light radius. Second, the trends of ordered to random motion ($(V/\sigma_{V,0})_{2.2}$) with stellar mass and sSFR (as a proxy for gas fraction) suggest that, on average, massive galaxies have marginally stable disks (Toomre, 1964) by $z\sim 2$. The higher gas fractions in galaxies at earlier epochs would naturally lead to more turbulent, thick disks under the Toomre disk stability criterion. This average gas fraction evolution may be enough to explain the observed structural evolution, where disks over similar parameter space are thicker at earlier times.

To disentangle the potential biases in this kinematics analysis and fully understand the source of the dynamical-baryonic mass tension for the higher-redshift galaxies, it is crucial to both directly measure molecular gas masses and have more detailed kinematic constraints. We will begin to address these open areas by examining molecular gas masses from the PHIBSS survey (Tacconi et al., 2018) for the MOSDEF galaxies that fall within both samples. Additionally, we will compare MOSDEF masses and kinematics with measurements of the same galaxies from the KMOS${}^{\hbox{\scriptsize{3D}}}$ survey (Wisnioski et al., 2015), to understand the limitations of the slit kinematic observations. Beyond these comparisons, future deep and high spatial resolution kinematic observations are also key to fully characterize the internal structures of high-redshift galaxies, particularly for small or low-mass galaxies, where the current kinematic constraints are most uncertain.

6 Conclusions

In this paper, we use spectra from the MOSDEF survey together with CANDELS HST/F160W imaging to study the kinematics and dynamical masses of 681 galaxies at $1.34\leq z\leq 3.8$, with stellar masses ranging from $M_{*}\sim 10^{9}M_{\odot}$ to $M_{*}\sim 10^{11.5}M_{\odot}$. In addition to kinematics and structural parameters, we use stellar masses derived from multiwavelength photometry and infer gas masses from either dust-corrected $\textrm{H}\alpha$ or $\textrm{H}\beta$ SFRs or SED SFRs if Balmer lines are unavailable.

We use the 3D kinematic models (MISFIT) developed in Price et al. (2016) to measure the galaxy kinematics from the misaligned galaxy-slit MOSFIRE observations. We use these models to measure both rotation and velocity dispersions from the 2D spectra for the 105 galaxies that have robust rotation detections. For the remaining 576 galaxies, we measure galaxy-integrated 1D velocity dispersions and use the kinematic models to convert the observed velocity dispersions into combined kinematic rms velocities and dynamical masses, assuming a fixed ratio $(V/\sigma_{V,0})_{R_{E}}$. These $(V/\sigma_{V,0})_{R_{E}}$ are derived as ensemble averages for the galaxies without detected rotation within bins of redshift, stellar mass, sSFR, $R_{E}$, and $\Sigma_{\mathrm{bar},e}$ by removing any trend between $\log_{10}(M_{\mathrm{dyn}}/M_{\mathrm{baryon}})$ and axis ratio $b/a$ (i.e., inclination).

We explore the relation between the ratio of rotation to velocity dispersion, $(V/\sigma_{V,0})_{2.2}$, as a function of stellar mass, sSFR, $\Sigma_{\mathrm{bar},e}$, and redshift. We find that $(V/\sigma_{V,0})_{2.2}$ increases with increasing stellar mass and decreases with increasing sSFR. These trends may indicate that these galaxies are marginally stable (Toomre, 1964), where galaxies with lower gas fractions (e.g., lower sSFR or higher stellar mass) will naturally have less turbulent motions. We additionally examine the relation between $V_{2.2}$ and $\sigma_{V,0}$ and stellar mass, sSFR, $\Sigma_{\mathrm{bar},e}$, and redshift for the galaxies with detected rotation. We find that $V_{2.2}$ is most correlated with $M_{*}$, and that $\sigma_{V,0}$ is most correlated with $\Sigma_{\mathrm{bar},e}$. There are weak, less significant trends between $\sigma_{V,0}$-$M_{*}$ and $V_{2.2}$-sSFR, and no correlation between $\sigma_{V,0}$ and sSFR. At fixed stellar mass, we see suggestive trends of $V_{2.2}$ increasing and $\sigma_{V,0}$ decreasing over time, but these trends are also generally consistent with no evolution given the uncertainties.

Using the mass measurements, we find that the median offset between dynamical and baryonic masses, $\Delta\log_{10}M=\log_{10}(M_{\mathrm{dyn}}/M_{\mathrm{baryon}})$, decreases with increasing redshift. The offset is relatively constant with stellar mass at fixed redshift. In contrast, we find that both larger galaxies and galaxies with lower surface densities tend to have higher mass offsets (i.e., lower $f_{\mathrm{bar}}$). The observed mass offset evolution implies an evolving dark matter fraction, where galaxies at $z\gtrsim 2$ are very strongly baryon dominated within their effective radii. The evolution of $f_{\mathrm{bar}}$ does not appear to be controlled by a single galaxy property, but instead reflects the intertwined effects of galaxy mass-size growth, gas fraction, and halo growth and evolution.

However, we find tension between the dynamical and baryonic masses at $z\gtrsim 2$, particularly for galaxies with small sizes or high densities. For these galaxies, the measured baryonic masses exceed the estimated dynamical masses (i.e., $M_{\mathrm{dyn}}<M_{\mathrm{baryon}}$). This mass discrepancy could be explained for a number of these galaxies having an offset Schmidt relation with higher SFEs (i.e., a smaller $\Sigma_{\mathrm{gas},e}$ can sustain the same $\Sigma_{\mathrm{SFR}}$), or if the galaxies have a higher virial coefficient $k_{\mathrm{tot}}$ (i.e., more spherical mass distributions). Nonetheless, our conclusions that galaxies become more baryon dominated toward higher redshifts and that massive galaxies generally have marginally stable disks by $z\sim 2$ are fairly robust even in the face of these measurement tensions.

The approach of using multiplexing, seeing-limited NIR spectrographs to constrain the average properties of high-redshift galaxies allows us to study kinematics for one of the largest samples of star-forming galaxies at $z\sim 1.5-3$, using a homogenous data set to extend down to much lower stellar masses than other surveys over a large range in redshift. However, further observations are necessary to reconcile the tension between the dynamical and baryonic masses at high redshifts and small sizes. In particular, direct observations of molecular gas masses are necessary to accurately measure baryonic masses. Detailed follow-up observations with adaptive-optics-assisted integral field unit (IFU) spectrographs are also crucial to better constrain the dynamical structures of these galaxies with unphysical baryon fractions.

We acknowledge useful discussions with T. T. Shimizu, S. Wuyts, E. J. Nelson, R. Genzel, L. Tacconi, S. Belli, K. Suess, and R. Trainor. We are grateful to the anonymous referee for their valuable comments. We thank the 3D-HST Collaboration, which provided the spectroscopic and photometric catalogs used to select the MOSDEF targets and derive stellar population parameters. This paper is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE 1106400 (S.H.P.). We acknowledge support from NSF AAG grants AST-1312780, 1312547, 1312764, and 1313171; archival grant AR-13907 provided by NASA through the Space Telescope Science Institute; and grant NNX16AF54G from the NASA ADAP program. The data presented in this paper were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. We extend special thanks to those of Hawaiian ancestry on whose sacred mountain we are privileged to be guests. We are most fortunate to have the opportunity to conduct observations from this mountain. This work is also based on observations made with the NASA/ESA Hubble Space Telescope, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. Observations associated with the following GO and GTO programs were used: 12063, 12440, 12442, 12443, 12444, 12445, 12060, 12061, 12062, 12064 (PI: Faber); 12177 and 12328 (PI: van Dokkum); 12461 and 12099 (PI: Riess); 11600 (PI: Weiner); 9425 and 9583 (PI: Giavalisco); 12190 (PI: Koekemoer); 11359 and 11360 (PI: O’Connell); 11563 (PI: Illingworth). This research made use of Astropy,666http://www.astropy.org a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013, 2018).  

Appendix A Kinematic Recovery Tests with MISFIT

Kinematics derived from slit spectra where the galaxy major axis and slit are misaligned suffer from a number of observational challenges. Any rotation signal will be blurred, as velocities from multiple radii are sampled within a given spatial slice, and some kinematic signal will be lost as portions of the galaxy fall outside of the slit. We account for the effects of galaxy-slit misalignments by fitting the kinematics using the forward models of MISFIT, which we developed and first presented in Price et al. (2016) (see also Section 3.1). However, the combined effects of galaxy size, slit misalignment, inclination, PSF FWHM, and S/N complicate the recovery of the intrinsic kinematics and thus may introduce scatter or bias in population-wide trends. In this appendix, we fit mock galaxy spectra with MISFIT and examine how well kinematics are recovered under a variety of observational conditions, in order to quantify the scatter and bias in our kinematic measurements. Fixed kinematic and structural parameters (with the exception of $R_{E}$) are selected for this test, to focus on the impact of the observational parameters (and to avoid prohibitively large mock data sets with long computation times).

We first generate a set of 1000 mock galaxies at $z=2$, each with $V_{a}=100\mathrm{\,km/s}$, $\sigma_{V,0}=50\mathrm{\,km/s}$, turnover radius $r_{t}=0.\arcsec 2$, Sérsic index $n=1.8$ (typical for emission line profiles at $z\sim 0.7-1.5$; Nelson et al. 2016), and intrinsic axis ratio $q_{0}=0.19$. We assume an instrument resolution of $\sigma_{\mathrm{inst}}=40\mathrm{\,km/s}$, which is typical for MOSFIRE in the K band (observing $\textrm{H}\alpha$ at $z\sim 2$). The following parameters are then randomly drawn from uniform ranges: galaxy-slit misalignment $\Delta\mathrm{PA}\in[-90^{\circ},90^{\circ}]$, inclination $i\in[0^{\circ},90^{\circ}]$ (face-on to edge-on), effective radius $R_{E}\in[0.\arcsec 1,1.\arcsec 25]$, PSF FWHM $\in[0.\arcsec 4,1.\arcsec 0]$, and central S/N per pixel $\in[3,20]$. The kinematic parameters are chosen to be fairly typical for our sample. The free parameter intervals (S/N, $\mathrm{FWHM_{PSF}}$, $R_{E}$, $\Delta\mathrm{PA}$, and $i$) are selected to cover the general range of values seen in the MOSDEF sample.

For each mock galaxy, we construct a model spectrum using the given fixed and randomly drawn structural and kinematic parameters, assuming that the PSF is Gaussian. Random noise is then added so that the central S/N per pixel equals the randomly selected S/N. Next, we use MISFIT to analyze the kinematics for each mock realization, using the same procedure used in the main analysis (see Section 3.1). From these results, we derive the best-fit $V_{2.2}$, $\sigma_{V,0}$, and $(V/\sigma_{V,0})_{2.2}$, and calculate $M_{\mathrm{dyn}}$ (using $V(R_{E})$, $\sigma_{V,0}$, and $R_{E}$) as in Section 3.3. Each mock observation is then classified as spatially resolved/aligned with detected rotation or with rotation limits, or as unresolved/misaligned following the criteria given in Section 3.1.

We show the offset between the input and recovered values of $V_{2.2}$, $\sigma_{V,0}$, $(V/\sigma_{V,0})_{2.2}$, and $M_{\mathrm{dyn}}$ (top to bottom) versus central S/N per pixel, PSF FWHM, emission line resolution $D_{\mathrm{50,emis}}/\mathrm{FWHM_{PSF}}$, galaxy-slit misalignment $|\Delta\mathrm{PA}|$, and inclination $i$ (left to right) in Figure 10. Overall, the parameters are recovered well on average, with little bias between the input and recovered parameters for both the full sample and for the subsample of resolved/aligned and rotation-detected mock observations. While the values are recovered well on average, there are relatively large scatters in the recovered $V_{2.2}$, $(V/\sigma_{V,0})_{2.2}$, and $M_{\mathrm{dyn}}$ for the full mock set ($\sim 0.25-0.4\mathrm{\,dex}$), though there is a lower scatter in the recovered $\sigma_{V,0}$ ($0.07\mathrm{\,dex}$). The scatter in the recovered parameters is lower when considering only the subset of mock observations that meet the resolution/alignment and detection cuts ($\sim 0.1-0.15\mathrm{\,dex}$ for $V_{2.2}$, $(V/\sigma_{V,0})_{2.2}$, and $M_{\mathrm{dyn}}$, and $0.04\mathrm{\,dex}$ for $\sigma_{V,0}$).

We see no systematic trends in the parameter recovery with central S/N per pixel or PSF FWHM. We see a slight trend of lower $V_{2.2}$ and $(V/\sigma_{V,0})_{2.2}$ when the emission line is barely resolved or unresolved (i.e., low $D_{\mathrm{50,emis}}/\mathrm{FWHM_{PSF}}$), so we thus restrict 2D kinematic fitting to galaxies with $D_{\mathrm{50,emis}}/\mathrm{FWHM_{PSF}}\geq 0.8$ (similar to the recovery results of Simons et al. 2016). There is a noticeable trend where $V_{2.2}$, $(V/\sigma_{V,0})_{2.2}$, and $M_{\mathrm{dyn}}$ are overestimated for large slit-galaxy misalignments, justifying the alignment cut where we do not fit 2D kinematics when $|\Delta\mathrm{PA}|>45^{\circ}$. We note that we still see a slight trend of lower $V_{2.2}$ and $(V/\sigma_{V,0})_{2.2}$ for the objects in the resolved/aligned and rotation-detected subsample with the smallest $R_{E}$ (maximum of $\sim\!0.1\mathrm{\,dex}$), but that there are no trends in the recovered $M_{\mathrm{dyn}}$ or $\sigma_{V,0}$. Finally, we see slight overestimates of $V_{2.2}$, $(V/\sigma_{V,0})_{2.2}$, and $M_{\mathrm{dyn}}$ for more face-on mock objects. However, this trend is not seen for the resolved/aligned and rotation-detected sample, as the kinematic quality cut appears to remove the objects most impacted by the face-on uncertainties.

Overall, this test demonstrates that while kinematics of individual objects may not be precisely recovered (resulting in scatter that may impact the recovery of trends when using individual measurements), the ensemble sample kinematic properties are well recovered with MISFIT, and thus our results do not suffer greatly from measurement biases.

Appendix B Comparison of Mixed SFR Indicators and Kinematics from Different Emission Lines

Successfully comparing trends of structure, kinematics, and matter content between different cosmic epochs requires either using the same measurement tracers or characterizing any bias between the tracers. In this appendix, we investigate whether biases arise in our analysis from using a combination of different SFR indicators and kinematics from different emission lines.

In this analysis, we use a “ladder” of SFR indicators, preferring Balmer-decrement-corrected $\textrm{H}\alpha$ SFRs when available (320 galaxies), but otherwise using (in order of preference) $A_{\mathrm{V}}$-corrected $\textrm{H}\alpha$ SFRs (176 galaxies), $A_{\mathrm{V}}$-corrected $\textrm{H}\beta$ SFRs (107 galaxies), or SED SFRs (78 galaxies), with $A_{\mathrm{V}}=A_{\mathrm{V,\,neb}}$ inferred from SED fit $A_{\mathrm{V,\,cont}}$ values. To test the accuracy of the lower-priority SFR indicators relative to the best-available Balmer-decrement-corrected $\textrm{H}\alpha$ SFRs, we select galaxies at $z\sim 2$ that have secure ($\mathrm{S/N}\geq 3$) $\textrm{H}\alpha$ and $\textrm{H}\beta$ detections. We restrict this test to only galaxies at $z\sim 2$ to provide the closest-available analogs to galaxies at $z\sim 3$, where only $A_{\mathrm{V}}$-corrected $\textrm{H}\beta$ and SED SFRs are available. We then alternatively derive the $A_{\mathrm{V}}$-corrected $\textrm{H}\alpha$ and $\textrm{H}\beta$ SFRs, as well as the SED SFRs, using the FAST SED fit results (see Section 2.1). The offset between each of these lower-priority SFR indicators and the Balmer-decrement-corrected $\textrm{H}\alpha$ SFRs is shown in Figure 11.

We find excellent agreement between the $A_{\mathrm{V}}$-corrected $\textrm{H}\alpha$ and $\textrm{H}\beta$ SFRs and the Balmer-decrement-corrected $\textrm{H}\alpha$ SFRs ($\Delta\log_{10}\mathrm{SFR}\sim 0.01-0.03\mathrm{\,dex}$), though there is a fair amount of scatter between the indicators ($\sigma_{{\mathrm{rms}}}\sim 0.2-0.3\mathrm{\,dex}$). However, we find that the SED SFRs are lower than the Balmer-decrement-corrected $\textrm{H}\alpha$ SFRs by $-0.23\mathrm{\,dex}$, with a moderately large scatter of $0.34\mathrm{\,dex}$. This offset implies that the SED SFR-inferred $M_{\mathrm{gas}}$ and $M_{\mathrm{baryon}}$ would underestimate the values derived from Balmer-decrement-corrected $\textrm{H}\alpha$ SFRs. Thus, accounting for this offset would lead to higher $M_{\mathrm{baryon}}$, increasing the tension between $M_{\mathrm{baryon}}$ and $M_{\mathrm{dyn}}$. Therefore, this bias cannot explain the observed tension between the dynamical and baryonic masses, but would rather amplify the problem. We do note that intrinsic differences between galaxies at $z\sim 2$ and $z\sim 3$ (e.g., metallicity, star-formation history shape), as well as differences in the available rest-frame photometry, may impact the accuracy of SED modeling, which could potentially lead to cases at $z\sim 3$ where SED SFRs actually overestimate $\textrm{H}\alpha$ SFRs (i.e., the very high sSFR objects at $z\sim 3$; Figure 1). If objects had higher SED than $\textrm{H}\alpha$ SFRs, $M_{\mathrm{baryon}}$ and the baryon fraction would be overestimated when using SED SFRs, which could potentially explain some of the dynamical-baryon mass tension. Also, the scatter between these indicators could potentially affect the quantification of trends involving individual galaxies and not just median values.

Next, we turn to potential biases from kinematic measurements with different emission lines. In our analysis, we use a mix of 2D and 1D measurements from $\textrm{H}\alpha$ (481 galaxies), [Oiii] (195 galaxies), and $\textrm{H}\beta$ (5 galaxies). Thus, we select objects for which we have measurements of both [Oiii] and $\textrm{H}\alpha$ or $\textrm{H}\beta$ and $\textrm{H}\alpha$. As with the SFR method comparison, we select galaxies at $z\sim 2$ to provide the best analogs to the $z\sim 3$ subsample. For this test, we do not require all three lines to be detected for each object, as the spectral coverage and skyline contamination of our spectra make clean kinematic measurements of all three lines in a single object relatively rare.

We first compare the 1D kinematics from different lines, as shown in Figure 12. The $\sigma_{V,\,\mathrm{obs}}$ values from $\textrm{H}\alpha$, [Oiii], and $\textrm{H}\beta$ agree very well on average, with offsets of only $-0.02\mathrm{\,dex}$ between [Oiii] or $\textrm{H}\beta$ and $\textrm{H}\alpha$. We do find some scatter ($\sigma_{\mathrm{rms}}\!\sim 0.1\mathrm{\,dex}$ and $\sim 0.2\mathrm{\,dex}$ for [Oiii] vs. $\textrm{H}\alpha$ and $\textrm{H}\beta$ vs. $\textrm{H}\alpha$, respectively), but overall mixing 1D dispersions measured from different emission lines likely does not introduce any bias into our analysis.

Finally, we compare the 2D kinematics measured from $\textrm{H}\alpha$, [Oiii], and $\textrm{H}\beta$. In Figure 13, we compare the values of $V_{2.2}$, $\sigma_{V,0}$, and $M_{\mathrm{dyn}}$ measured from 2D fitting of [Oiii] or $\textrm{H}\beta$ with those from $\textrm{H}\alpha$. The values measured from $\textrm{H}\alpha$ and the other lines are generally in good agreement, with at most a $\sim 0.13\mathrm{\,dex}$ offset. We do find a fair amount of scatter in the recovered values, but this is generally smaller than the uncertainty in the offsets. Nonetheless, as with the kinematic recovery and SFR indicator tests, the scatter between the different emission lines’ 2D kinematics could possibly influence the investigated galaxy property trends. Although the 2D comparison sample is small, this test suggests that the results derived from different emission lines are compatible, so that mixing of kinematic tracers does not systematically impact our results.